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Motivation 


This  work  addresses  the  needs  for  advanced  computational  tools  and  theoretical 
models  devoted  to  the  study  of  two  interrelated  topics  of  interest  to  the  Air  Force: 
Radiation  Belt  Remediation  (RBR)  and  potential  use  of  tethers  as  low  frequency 
oscillators/radiators.  Both  problems  stretch  the  capabilities  of  existing  numerical 
methods,  while  connecting  with  a  body  of  previous  work  on  probe  theory,  spacecraft- 
plasma  interactions  and  electrodynamic  tether  propulsion.  The  Electrostatic  RBR 
application  would  typically  require  MV-level  potentials  and  complex,  multi-strand  wire 
arrangements.  The  high  voltages  imply  extreme  disparity  of  length  scales,  as  well  as 
relativistic  conditions  in  some  cases.  The  geometrical  complexity  and  driving  physical 
processes  may  require  3D  capabilities,  particularly  when  magnetic  effects  cannot  be 
ignored.  For  its  part,  the  study  of  sheath  ion  or  electron  oscillations  in  the  vicinity  of  a 
high  power  radiating  tether  requires  tracking  of  a  time -dependent  sheath  boundary  and 
use  of  boundary  conditions  that  allow  radiation  escape  while  denying  spurious 
reflections. 

Our  studies  also  addressed  the  engineering  integration  aspects  of  these  two  problems. 
The  detailed  kinetic  code,  which  is  required  for  accurate  and  consistent  modeling  of  the 
physics  is  expected  to  be  ill-suited  for  parametric  studies  that  aim  at  design  and 
optimization  of  an  actual  RBR  or  broadcasting  system.  For  this  reason,  we  also 
developed  reduced  order  models  that  retain  the  main  scaling  and  performance  capabilities 
while  drastically  reducing  computational  effort.  The  applications  are  in  principle 
physically  different,  so  we  developed  separate  models  for  the  electrostatic  and  for  the 
wave-driven  cases.  These  models  can  be  ultimately  used  to  guide  the  selection  of 
parameters  for  detailed  simulations,  as  well  as  for  a  future  systems  study  that  will 
produce  recommendations  regarding  system  configuration  and  project  power,  mass  and 
other  global  metrics  for  each  application. 


Background 

The  RBR  problem  arises  from  the  possibility  of  man-made  intensification  of  the 
already  dangerous  levels  of  high-energy  particle  populations  trapped  by  the  geomagnetic 
field  in  the  altitude  range  of  1,000-15,000  km.  These  particles,  typically  protons  and 
electrons,  have  pitch  angles  with  respect  to  the  magnetic  field  high  enough  to  induce 
repeated  reflections  above  the  dense  atmosphere  altitude.  Any  intervention  that  can 
reduce  this  pitch  angle  can  ultimately  lead  to  absorption  of  the  particle  as  it  dips  farther 
into  the  atmosphere  near  magnetic  poles.  Among  the  possible  mechanisms  that  have  been 
proposed  for  this  purpose,  two  deserve  mention  here: 

(a)  injection  of  whistler  waves  into  the  magnetospheric  cavity  to  mimic  the  particle¬ 
scattering  action  of  ionospheric  hiss  or  lightning,  which  are  some  of  the  natural  decay 
mechanisms  for  the  trapped  particles,  and 

(b)  scattering  by  some  high  potential  structure  (a  metallic  tether  or  combination  of 
tethers  has  been  suggested)  that  will  send  a  fraction  of  the  protons  into  the  loss  cone  for 
their  magnetic  line. 
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Our  study  addresses  both  of  these  mechanisms.  In  addition  to  the  direct  interest  for 
this  RBR  application,  our  studies  will  have  applicability  to  probe  theory,  electrodynamic 
tether  propulsion  and  power  generation,  and  spacecraft-plasma  interactions. 

The  study  of  RBR  tether  operation  can  be  attempted  at  two  levels  of  accuracy.  For 
design  studies,  where  multiple  parametric  trials  need  to  be  evaluated,  a  simple  model  can 
be  developed  that  uses  an  approximate  analytical  potential  distribution  around  the  tether, 
evaluates  the  scattering  cross-section  for  particle  removal,  and  uses  it  to  calculated  total 
loss  rate.  On  the  other  hand,  for  a  clear  understanding  of  he  physical  issues  and  for 
numerical  precision,  one  needs  to  resort  to  a  detailed  kinetic  simulation  capable  of 
including  magnetic  effects,  relativistic  corrections  and  complicated  geometries. 


For  the  scattering  approach,  a  good  example  of  the  simplified  approach  is  the  original 
study  of  Danilov  et  al  [1],  Their  work  can  be  used  to  establish  some  of  the  main 
parameters  of  a  typical  RBR  system,  and  this  can  then  guide  our  discussion  of  the 
numerical  issues  for  the  kinetic  simulation.  Danilov  pointed  out  that  a  pair  of  tethers 
biased  in  a  double-probe  arrangement  could  be  used,  with  most  of  the  potential  deviation 
from  the  background  being  taken  up  by  the  negative  tether,  which  then  becomes  the 
active  scattering  center.  Since  the  trapped  proton  energy  is  up  to  1  MeV,  the  bias  must 
also  be  of  this  order.  The  electron  temperature  in  the  background  magnetospheric  plasma 
is  of  the  order  of  100  eV,  and  the  plasma  density  at  the  high  altitudes  considered  can  be 
as  low  as  10A8  m-3,  for  a  Debye  length  of  nearly  25  m.  The  radius  of  the  sheath  formed 
around  the  negative  tether  can  be  estimated  [Ref  2]  from  the  implicit  formula 


p  v  v 

—  =  2.554(-^> 
kT,  dD 


(1) 


where  V  is  the  tether  potential,  Te  the  electron  temperature,  rs  the  sheath  radius,  R  the 
tether  radius,  and  do  the  Debye  length.  For  the  conditions  quoted  above,  with  a  1MV 
potential,  the  sheath  radius  ranges  from  1,400  m  when  the  tether  radius  is  0.01  mm  to 
1,900  m  when  the  tether  radius  is  10  mm.  The  ratio  of  sheath  to  tether  radii  is  then  from 
1.4xl08  to  1.9xl05.  The  energetic  particle  density  is  much  smaller  than  the  background 
density,  which  sets  the  Debye  length  and,  together  with  the  tether  radius,  determines  the 
current  capture  rate  by  the  two  tethers,  and  hence  the  power  required  by  the  system. 


Numerical  work  towards  the  RBR  problem  has  been  carried  out  by  Choiniere  and 
Gilchrist  [2]  as  well  as  by  Minor,  of  Tethers  Unlimited  (Ref.[3]).  The  result  expressed  by 
Eq  (1)  is  a  correlation  of  the  numerical  results  of  Ref.  [2],  and  has  since  been  confirmed 
through  analytical  work  of  Sanmartin  et  al.  [4],  References  [2,3]  have  used  an 
axisymmetric  steady  state  model  in  which  the  potential  is  gradually  adjusted  as  a  result  of 
a  large  number  of  charged  particle  trajectories  that  are  computed  at  each  stage  in  the 
iteration.  This  procedure,  analogous  to  that  used  in  the  classical  probe  studies  of 
Laframboise  [5],  accelerates  convergence  at  the  cost  of  ignoring  the  transient  dynamics. 
In  addition  to  its  direct  importance  for  oscillator  studies,  this  transient  behavior  may  also 
be  relevant  to  the  scattering  problem  if,  as  suggested  by  our  electron  capture  result  [6],  it 
results  in  the  long-term  trapping  of  a  significant  population  around  the  attracting  tether. 
The  use  of  an  electrostatic  tether  as  a  combination  resonator-antenna  for  broadcasting  of 
VLF  or  ELF  signals  appears  to  have  received  comparatively  less  attention  than  the  RBR 
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application,  although  Refs.  [17-19]  have  studied  the  physical  basis  and  considered  the 
injection  of  whistlers  from  antennas  in  orbiting  spacecraft.  Scaling  arguments  were 
advanced  by  NRL  researchers  (Ref.  [7]),  who  pointed  out  that  ion  oscillations  in  the 
sheath  of  a  negative  tether  would  have  frequencies  in  the  range  of  interest.  Whether  these 
oscillations  can  become  self-sustained,  or,  on  the  contrary,  will  require  steady  excitation 
power,  is  at  this  time  unclear.  Transients  in  the  TSS-1  space  experiment  were  described 
by  Bilen  et  al.  (Ref.  [8])  and  an  experimental  and  numerical  study  of  the  possibility  of 
enhancing  electron  current  collection  by  the  scattering  effect  of  forced  sheath  potential 
fluctuations  was  reported  by  Choiniere  et  al.  (Ref.  [9]).  The  effect  was  indeed  positive, 
especially  near  resonance  (somewhat  below  plasma  frequency),  but  the  driving  power 
was  also  very  high  at  that  frequency.  Current  enhancement  per  unit  power  peaked  well 
below  resonance. 

One  of  the  most  likely  sources  of  oscillation  is  the  trapped  population  of  ions  that  can 
develop  around  a  negatively  biased  tether.  The  possibility  of  these  trapped  populations  in 
potential  wells  has  been  long  recognized  (Ref  [10]).  Gurevich  [11]  analyzed  their 
development  though  adiabatic  potential  changes  in  one-dimensional  geometries,  and 
derived  the  resulting  energy  distribution  of  the  trapped  particles.  He  also  recognized  that 
fast  transients  could  result  in  very  similar  distributions  as  well.  Sanmartin  et  al  [12] 
generalized  Gurevich’s  analysis  to  axisymmetric  geometries  and  explained  through  the 
presence  of  trapped  electrons  the  paradox  that  arises  in  the  frontal  pre-sheath  region  of  a 
positively  biased  tether  in  a  hypersonic  flow:  the  retarded  ions  tend  to  pile  up  and 
overshoot  the  upstream  plasma  density,  while  the  untrapped  electron  density  cannot 
exceed  that  upstream  density  [13],  In  our  electron  capture  PIC  simulations  [5],  we  have 
observed  large-scale  trapping  of  electrons  inside  the  sheath  around  a  positive  tether,  as 
well  as  weaker  trapping  in  the  pre-sheath.  Their  origin  is  clearly  the  potential  transient  at 
tether  turn-on,  but  we  have  also  noticed  that  there  is  a  sort  of  electrostatic  feedback 
tending  to  restore  that  population  when  scattering  or  artificial  removal  is  used  in  the 
simulation  to  reduce  it. 


Given  the  presence  of  a  charge  population  orbiting  a  biased  tether,  it  is  not  difficult  to 
conceive  the  possibility  of  “sloshing”  modes  of  this  charge,  and  we  have  occasionally 
noticed  these  modes  in  simulations  (probably  triggered  by  numerical  events).  Their 
frequency  is  of  the  order  of,  but  lower  than  the  plasma  frequency  of  the  trapped  type  of 
particles.  This  is  supported  by  the  following  simplified  argument:  it  can  be  easily  shown 
that  the  frequency  of  a  harmonic  oscillator  in  which  the  peak  potential  energy  Umax 

I  ZU 

occurs  at  maximum  deflection  xmax  is  given  by  co  =  J — .  For  the  particles  trapped 

V  MX  max 

inside  the  sheath,  the  potential  energy  is  a  sizeable  fraction  of  eV,  while  the  amplitude  of 
the  oscillation  can  be  a  similar  fraction  of  rs,  the  sheath  radius.  Equation  (1)  shows  an 
approximate  scaling  of  rs  as  (eV/kTe)3/4.  Combining  these  expressions,  one  can  show  that 


% 

CO  « - - — 7-J 

0 eV/kTe)u 4 


(2) 


where  oop  is  the  plasma  frequency  of  the  oscillating  particle  (ions  in  the  case  of  present 
interest). 
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RBR  Research  Requirements 

The  discussion  in  previous  section  can  serve  as  a  basis  for  identification  of  the 
research  needs  and  challenges  in  this  area.  Some  of  these  have  been  addressed  and 
resolved  in  our  work,  while  others  remain  as  challenges  for  future  research.  Regarding 
the  high-voltage  tether  option,  the  very  large  sheath  radius  makes  it  difficult  to  devise  a 
numerical  scheme  that  can  resolve  both,  the  thin  tether  neighborhood,  and  the  sheath 
edge,  which  because  of  its  large  local  potential  gradients,  may  significantly  influence 
scattering.  A  possible  compromise  would  involve  the  use  of  an  “analytical  inner  domain”, 
as  in  the  work  of  Onishi  [14],  dynamically  coupled  to  the  coarser  outer  grid.  The  inner 
domain  is  essentially  Coulombic  and  axisymmetric,  while  the  outer  domain  needs  to 
account  for  the  relative  plasma  flow  and  the  geomagnetic  field,  and  has  to  be  at  least  two- 
dimensional,  with  3D  vector  quantities.  Since  many  of  the  promising  designs  include 
multiple  tethers,  these  outer  domains  need  to  include  more  than  one  sheath,  and  the 
numerical  domain  will  have  to  extend  for  several  hundred  Debye  lengths.  Because  the 
spatial  grid  for  solution  of  the  Poisson  equation  needs  to  resolve  the  Debye  length,  one 
further  difficulty  accrues  in  that  a  very  large  number  of  grid  nodes  will  be  necessary  even 
for  these  outer  regions  alone.  A  possible  approach  that  can  be  envisioned  is  a  transition  to 
a  quasi-neutral  fluid  model  of  the  electrons  at  some  distance  from  each  tether;  the  Poisson 
equation  then  degenerates  to  quasineutrality,  ne=rii,  and  the  grid  restriction  to  Debye  size 
disappears.  We  have  experience  with  this  approach  in  the  plasma  plume  area  (Ref.  [15]). 

For  the  oscillator  problem,  one  faces  the  need  to  resolve  a  shifting  sheath  boundary, 
probably  requiring  the  use  of  an  adaptive  grid  scheme.  Also,  since  the  emphasis  is  now 
on  a  consistent  rendering  of  the  oscillatory  motions,  the  boundary  conditions  used  at  the 
edges  of  the  computational  domain  will  need  to  be  of  the  radiation  type,  allowing 
outgoing  waves,  but  not  incoming  waves  to  cross.  This  is  a  familiar  problem  in  CFD  and 
acoustics  computations,  where  only  one  simple  type  of  wave  exists.  For  magnetized 
plasma,  care  will  have  to  be  exerted  to  accommodate  all  the  wave  modes  of  interest.  Our 
analytical  work  on  Whistler  emission  (Appendix  4)  can  be  used  for  this  task.  One 
additional  requirement  for  this  type  of  computation  is  the  accurate  determination  of  the 
trapped  particle  population,  for  which  not  much  numerical  experience  exists  as  yet.  At  a 
minimum,  the  numerical  scheme  will  need  to  conserve  energy  over  the  full  computational 
time,  rather  than  just  the  passage  time  of  a  non-trapped  particle,  as  is  more  usually  the 
case.  Perhaps  the  feedback  mechanism  referred  to  above,  if  fully  confirmed,  can  be  relied 
on  to  preserve  the  trapped  population  even  in  the  presence  of  long-term  numerical  error 
accumulation 
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Reported  Results 

During  the  period  of  this  Grant  work  was  done  along  several  paths:  collection  and 
review  of  experimental  data,  reduced  order  model  studies  of  the  electrostatic  RBR 
problem,  theoretical  work  on  remediation  using  plasma  waves,  and  development  of 
relativistic  adaptive  numerical  model  for  MV-tethers,  ambient  plasma  dynamics  with 
self-consistent  description  of  wave-plasma  interactions,  and  analytical  modeling  of  the 
radiation  pattern  from  a  long  tether  acting  as  a  VLF  antenna.  This  work  is  described 
below  in  summary  form,  and  additional  details  are  available  in  the  Appendices. 


Task  1.  Environment  Parameters  for  RBR 

Initial  analysis  included  review  of  Radiation  Belts  composition  data  from  20+  US  and 
Russian  sources,  which  are  partially  included  in  the  reference  list  at  the  bottom  of  this 
document,  with  additional  details  in  Appendix  1.  Charge  and  neutral  particle 
concentrations  and  energy  spectra  have  been  identified  for  various  spatial-temporal 
conditions,  etc. 

Next  a  comparison  of  natural  Van- Allen  belts  conditions  with  artificial  ionic  and 
electron  belts  was  performed  to  better  characterize  the  space  environment  to  be 
remediated  with  high-voltage  tethers  and  RF -waves.  Characterization  of  plasma 
collisionality,  Debye  shielding,  thermal  conditions  has  been  performed. 

Finally,  the  analysis  of  bare  tethers  erosion  loses  in  orbit  was  done  using  available 
sputtering  data.  In  short,  the  conclusion  is  that  this  effect  can  be  neglected  for  the  tether’s 
operational  parameters  currently  envisioned  for  electrostatic  RBR. 


Task  2.  Study  of  electrostatic  scattering  by  a  tether  as  an  RBR  system 

The  systems  model  for  radiation  belt  remediation  is  the  development  of  work  we  had 
initiated  prior  to  this  work  [16].  The  model  assumes  a  purely  Coulombic  potential  around 
a  negative  tether,  extending  to  zero  value  at  the  radius  rs  of  the  sheath.  The  tether  is 
assumed  to  be  in  a  circular  equatorial  orbit,  directed  vertically  (see  Fig.l).  Magnetic 
effects  are  ignored  during  the  scattering  interaction,  but  a  background  magnetic  field  is 
assumed,  with  its  dipole  axis  coincident  with  the  Geographic  NS  axis. 

A  classical  scattering  formulation  furnishes  a  relationship  between  the  incoming 
particle  “miss  distance”  in  the  equatorial  plane  and  its  angular  deflection  by  the 
electrostatic  potential  field.  Analytical  conditions  are  derived  to  test  whether  the  leaving 
direction  falls  within  the  Loss  Cone  at  the  location  of  the  interaction  (defined  as  the 
angular  range  around  the  magnetic  direction  that  will  allow  penetration  of  the  particle 
below  1 10  km  in  either  polar  region),  and  an  integration  is  performed  over  all  velocities 
belonging  to  a  “hollow  cone”  Maxwellian  distribution  and  over  all  miss  distances 
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consistent  with  entering  the  Loss  Cone.  The  calculations  are  only  partially  done 
analytically,  the  rest  being  performed  as  a  series  of  nested  numerical  quadratures. 


Figure  1  -  Geometry  of  RBR  setup 
used  for  reduced  order  model  analysis 


The  program  embodying  the  above  formulation  has  been  coded  and  tested. 
Extensions  are  added  to  calculate  the  current  collected  by  either  polarity  tether,  hence  the 
DC  power  required  to  sustain  the  potential.  The  OML  formulation  is  used  for  this 
purpose.  This  tool  was  exercised  to  postulate  and  evaluate  a  number  of  system 
configurations  (orbital  parameters,  number  of  tethers,  tether  dimensions,  system  mass, 
expected  life,  cost,  etc).  These  are  applied  to  assess  against  the  system’s  performance, 
namely,  the  time  to  reduce  the  initial  RB  energetic  particle  population  by  a  desired 
percentage.  Natural  as  well  as  enhanced  initial  conditions  are  considered  in  the  model, 
and  the  competition  of  the  natural  capture  rate  of  cosmic  radiation  is  an  input  to  the 
calculations.  A  detailed  description  of  the  reduced  model  can  be  found  in  MS  thesis  by 
Chris  Zeineh  [16],  which  is  incorporated  here  as  Appendix  2.  A  useful  presentation 
version  of  this  work  is  also  given  in  Appendix  5. 


Task  3.  Study  of  the  Wave-Based  RBR 


3.1  General  discussion 

The  loss  mechanisms  for  the  natural  population  of  energetic  particles  in  the  Van 
Allen  belts  have  been  discussed  quantitatively  by  Abel  and  Thorpe  [17,18],  and  later 
examined  in  relation  with  the  Radiation  Belt  Remediation  concept  by  Inan,  Bell  and 
Bortnik  [19].  Abel  and  Thorpe  concluded  that,  except  for  the  lowest  altitudes,  where 
Coulomb  scattering  dominates,  natural  and  artificially  injected  waves  of  the  Whistler 
family  are  responsible  for  most  of  the  precipitation  of  trapped  particles.  The  mechanism 
by  which  a  fraction  of  the  trapped  particles  gradually  lose  pitch  angle  and  eventually 
enter  the  loss  cone  appears  to  be  the  accumulation  of  short  periods  of  resonant  interaction 
between  a  given  trapped  particle  and  the  ducted  radiation  that  propagates  along  the  same 
magnetic  tube.  This  mechanism  is  represented  as  a  diffusion  in  pitch  angle  space,  and  the 
authors  present  a  quasi-linear  formulation  to  calculate  the  corresponding  diffusivities  due 
to  the  various  resonant  orders,  including  O’th  order  (Landau  resonance),  first  order 
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(cyclotron  resonance)  and,  at  the  higher  energies,  second  and  higher  orders  as  well.  The 
frequency  range  that  dominates  precipitation  depends  on  particle  energy  and  altitude  ( L - 
shell);  in  particular,  for  the  inner  belt  (L< 2.6),  the  major  contributor  was  found  to  be  the 
near-monochromatic  radiation  from  a  small  number  of  identifiable  YLF  antennas 
operating  around  the  globe.  This  finding  lends  credence  to  the  possibility  of  a  deliberate 
intervention  that  could  strongly  reduce  the  mean  life  of  these  particles,  possibly  well 
below  the  roughly  100  days  that  apply  to  the  inner  radiation  peak. 

Inan,  Bell  and  Bortnik,  in  their  later  work  [19],  discussed  several  parametric  changes 
that  would  lead  to  further  acceleration  of  the  pitch  angle  diffusion.  Thus,  it  was  shown 
that  VLF/ELF  frequencies  of  only  a  few  kHz,  rather  than  the  roughly  20  kHz  accounted 
for  in  the  Abel  and  Thorpe  study,  would  accelerate  the  effect  by  factors  of  the  order  of 
30.  Additional  amplification  should  occur  through  the  efficient  accumulation  of  whistler 
wave  energy  in  the  magnetospheric  cavity. 

As  Inan  et  al.  point  out,  an  electrical  dipole  antenna  radiates  most  intensely  at  a 
frequency  for  which  the  wavelength  is  twice  the  antenna  length  ( L  «  c/(2nf)),  where  n  is 
the  index  of  refraction  at  a  frequency  /  and  at  the  particular  wave  direction  y/  with  respect 
to  the  basic  B  field.  Except  for  large  y/  angles,  near  the  resonance  angle  (ie,  near  90°),  n 
is  of  the  order  of  15,  so  that  a  choice  f=2.5  kHz,  as  proposed  by  Inan  el  al.  implies  an 
antenna  length  of  about  4  km.  It  seems  clear  that  conventional  satellite-based  antennae, 
limited  as  they  would  be  to  L  of  the  order  of  100m,  would  force  the  use  of  much  higher 
frequencies,  with  the  attendant  lower  scattering  effectiveness.  To  circumvent  this 
problem,  Inan  et  al.  propose  injection  angles  near  the  resonance  angle,  for  which  the 
index  of  refraction  can  be  as  high  as  1000  (limited  by  thermal  effects).  If  instead  one  used 
a  fraction  of  the  length  of  a  multi-km  gravity-gradient  stabilized  tether,  the  length  of  the 
antenna  could  be  easily  selected  to  match  any  desired  frequency,  down  to  below  1kHz,  as 
well  as  any  convenient  wave  angle.  Concepts  using  several  tethers  segments  as 
independent  dipoles  could  offer  additional  flexibility,  including  the  possibility  of  phasing 
arrangements  that  would  help  focus  the  radiated  power. 


3.2  Pitch  scattering  by  Whistlers 

We  have  pursued  modeling  work  along  the  lines  of  Ref.  [17]  in  order  to  create  a  tool 
for  rapid  evaluation  of  a  variety  of  radiating  tether  concepts.  Our  model  incorporates  the 
effect  of  the  n  ^  0  interaction  modes  and  their  bounce-averaged  contributions  to  the  pitch 
diffusivity  as  a  function  of  the  electron  pitch  and  energy.  The  preliminary  result,  see  Fig. 
2,  agrees  well  with  those  reported  in  Ref  [17].  A  more  detailed  account  of  this  work  is 
given  in  our  Appendix  3,  and  in  presentation  form  in  Appendix  5. 
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Pitch  Angle  Diffussion  Coefficient 


Figure  2  -  Calculated  pitch-angle  diffusion 
coefficient  for  the  quasi-linear  model 


3.3  Whistler  radiation  from  a  tether 

Future  work  on  the  spatial  integration  for  a  given  L-shell  requires  as  inputs  the 
angular  antenna  pattern  and  a  measure  of  the  local  wave  intensity,  such  as  the  wave  B 
field  amplitude.  This  was  initially  attempted  by  using  the  results  by  Wang  and  Bell 
[20,21]  for  cold  and  warm  magneto-plasmas,  but  it  was  thought  important  to  develop  our 
own  detailed  model,  still  within  the  cold  plasma  linear  approximation,  so  that 
generalizations  to  other  wave  bands  or  other  geometries  could  later  be  easily 
accommodated.  This  work  is  nearing  completion  now,  as  the  MS  Thesis  for  Mr.  Yu 
Takiguchi,  and  its  current  status  is  described  in  detail  in  Appendix  4  (and  in  presentation 
form  in  Appendix  5).  Refinements  will  be  later  introduced  based  on  our  numerical  code 
(see  the  following  sections).  The  Thesis  is  expected  by  September  2008,  and  it  will 
contain  a  complete  theoretical  analysis  of  the  radiation  pattern  as  a  function  of  antenna 
(tether)  size  and  orientation,  plus  a  full  discussion  of  the  electron  dynamics  as  driven  by 
the  wave,  for  use  in  further  pitch  diffusion  calculations.  Numerical  results  will  be  shown 
for  parameter  ranges  of  interest  to  RBR.  Beyond  this,  we  intend  to  apply  the  same 
methodology  to  the  lower  frequency  ion  band,  which  is  the  one  likely  to  be  of  use  for 
removal  of  protons  and  other  positive  ions  from  the  Belts.  In  addition,  our  numerical 
methods  will  be  used  to  include  nonlinear  near-tether  effects  and  complex  geometries. 
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Task  4.  Development  of  an  Adaptive  Numerical  Model  for  RBR 

The  several  orders  of  magnitude  difference  in  spatial  and  temporal  scales  in  the  RBR 
electrostatic  tether  problem,  and  the  smaller,  but  still  large  variation  of  scales  for  the  low 
frequency  driven  electromagnetic  tether,  makes  impossible  the  direct  modeling  of  such 
systems  using  traditional  kinetic  approaches,  e.g.  a  regular  particle-in-cell  method. 

Indeed,  for  the  plasma  of  interest  the  ratio  of  Debye  length  to  the  tether  diameter 
could  be  on  the  order  of  104  and  more.  Unlike  common  electrodynamic  tethers,  with 
sheaths  sizes  comparable  to  the  radius  of  the  wire  and  to  the  Debye  length,  the  mega- volt 
electrostatic  tethers  have  sheath  sizes  in  the  excess  of  hundreds  of  meters,  comprising 
hundreds  of  Debye  lengths.  At  the  same  time  accurate  calculation  of  small-scale 
phenomena  are  required,  for  instance: 

i)  current  collection  by  thin  tether  (exact  trajectory  intersection  with  surface); 

ii)  parameters  of  the  scattering  cross-sections,  affected  by  space  charge  near  wire. 

The  simulation  domain  to  capture  the  important  pre-sheath  region  has  to  be  at  least 
several  sheath  lengths,  10  if  possible.  Thus  we  easily  arrive  at  a  billion  required  cell 
grids,  making  simulations  impractical  if  a  mostly  uniform  grid  is  retained,  because  10- 
100  particles  are  required  per  cell  to  assure  ~10%  statistical  accuracy  of  calculating 
transient  regimes.  This  problem  can  be  alleviated  by  either  by: 

i)  analytical  integration  of  the  trajectories  near  a  high-voltage  tether,  as  was  first 

proposed  by  Onishi  [14],  or 

ii)  numerical  integration  of  trajectories  in  the  high  potential  gradient  region  at  a  much 

finer  scale  [25], 

Because  the  singular  regions  amount  to  about  1%  of  the  entire  simulation  domain, 
both  procedures  are  acceptable  from  the  computational  standpoint.  However,  the  purely 
analytical  approach  does  not  allow  potential  variations,  e.g.  due  to  trapped  particle  space- 
charge,  while  the  numerical  approach  may  require  tiny  sub-stepping  near  the  wire. 

A  significant  acceleration  of  computations  could  be  obtained  by  using  a  variation  of 
FFT  technique  that  allows  internal  boundaries.  It  was  implemented  in  for  the  2D  case  in 
[25],  and  used  for  detailed  studies  of  current  collection  by  tethers  in  different 
configurations  [27,28],  including  multiple  tether  interaction.  But  this  method,  though 
used  in  many  cases,  still  relies  on  the  fact  that  the  grid  is  uniform,  and  therefore,  has  a 
limited  potential  to  study  tethers  with  dimensions  much  smaller  than  Debye  length,  and 
with  sophisticated  shapes. 

For  practical  engineering  purposes,  an  additional  important  requirement  is  the  ability 
to  simulate  complex  geometries  of  the  tethers,  especially  when  we  are  talking  about 
predictive  modeling  of  tether  arrays,  systems  with  multiple  elements  as  VLF  antennas 
with  complex  cross-sections,  etc.  Partial  patches,  as  shown  in  Fig.  3,  can  be  used  in  a 
limited  number  of  cases,  because  in  the  general  case  these  cross-sections  are  not  square 
or  circular,  and  it  is  impossible  to  match  an  arbitrary  shape  with  a  combination  of  fixed 
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elements.  Moreover,  tethers  can  move  mutually, 
excluding  fixed  grids  as  candidates  for  the 
universal  approach. 


Thus,  despite  the  possibility  of  using  some  of 
the  partial  remedies  so  far  discussed,  it  becomes 
obvious  that  the  best  way  to  proceed  with  multi¬ 
scale  modeling  of  systems  with  complex 
geometries  is  the  use  of  non-uniform  grids. 

These  are  most  desirable  if  they  have  an 
additional  ability  of  conforming  automatically  to 
the  arbitrarily  shaped  collecting  surfaces.  One  of 
such  possibility  is  illustrated  in  Fig.  4,  where  a 
circular  tether  is  approximated  using  uniform  Figure  3  -  Square-to-round  patch  for 

structured  and  non-uniform  unstructured  grids.  uniform  grid  near  the  tether  [14] 
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Figure  4  -  Tether  shape  approximation  using  uniform  structured 
(left)  and  non-uniform  unstructured  grid  (right  panel) 


Modeling  on  non-uniform  grids 
allows  capturing  of  the  tethered  RBR 
and  VLF  problems  with  just  millions  of 
cells,  instead  of  billions.  Moreover,  the 
adaptive  approach  opens  the  possibility 
to  study  interference  of  several  tethers 
without  doubling  the  amount  of 
calculation.  Indeed,  a  non-uniform 
adaptive  mesh  is  capable  of  refining  the 
mesh  near  the  wire,  leaving  the 
majority  of  the  nodes  large.  As  a  result, 
the  total  number  of  nodes  is  increased 
by  a  few  percent.  The  computational 
grid  in  the  tethers’  vicinity  is  refined,  as 
could  be  seen  from  Fig.  5,  where  an 
example  of  a  two-tether  system  is 


tethers 

\ 

\ 


\m 

In 

In 

Figure  5  -  Example  of  an  adaptive  grid 
local  refinement  for  a  2-wire  system 


12 


presented.  It  is  clear  that  vast  exterior  region  of  the  simulated  domain  remains 
unmodified.  Doubling  resolution  just  near  a  single  pole  requires  adding  only  4  extra  cells. 
To  reduce  resolution  by  106  ( — 220)  only  160  extra  nodes  will  be  required. 

We  can  easily  match  any  spatial  dimension  through  the  exponential  subdivision  of 
cells:  just  14  divisions  are  required  to  reduce  cell  size  in  excess  of  4  orders  of  magnitude. 
Estimates  show  that  simultaneous  capturing  on,  say,  10x1 0km  domain,  and  a  tether  with 
1mm  radius  will  require  about  0.5  million  cells,  and  about  20  millions  of  particles.  Such 
simulation  could  be  easily  carried  out  on  a  desktop  PC. 


The  adaptive  grid  approach  allows  not  only  flexible  internal,  but  also  external 
boundaries.  Usually  calculations  are  performed  in  a  rectangular  domain,  which  is  not 
always  the  best  choice.  For  instance,  a  circular  domain  allows  cutting  the  number  of 
required  cells  by  a  quarter,  as  could  be  seen  from  Fig.  7,  in  which  contours  of 
electrostatic  potential  are  presented  for  regions  with  two  different  shapes.  Also  it  shows 
initial  grids,  and  meshes,  after  they  automatically  adapted  to  the  structure  of  the  potential. 
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Figure  7  -  Comparison  of  two  electrostatic  tether 
problem  in  a  rectangular  and  in  a  circular  domain 

We  also  consider  the  introduction  of  other  computational  techniques  to  speed-up 
calculations.  As  an  important  example,  since  most  of  the  computational  cells  may  be 
located  in  smoothly  varying  distant  pre-sheath  regions,  we  will  evaluate  the  use  of  hybrid 
methods,  where  a  fluid  description  for  electrons  is  applied  for  these  quasi-neutral  parts  of 
the  domain.  These  regions  of  hybrid  simulation  will  dynamically  be  matched  to  those 
where  the  full  Poisson  equation  needs  to  be  solved,  as  we  did  in  Ref.  [15].  The  boundary 
is  adjusted  by  evaluating  the  degree  on  non-neutrality  implied  by  the  quasi-neutral 
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solution  in  its  vicinity.  The  potential  derived  from  the  quasi-neutral  side  becomes  a 
Dirichlet  condition  for  the  Poisson  side,  and  an  electron  injection  scheme  has  to  be 
implemented  on  the  kinetic  side  of  the  boundary.  The  results  of  these  hybrid  methods  are 
being  verified  against  fully  kinetic  modeling  for  selected  test  cases.  The  development  of 
the  latter  becomes  the  top  priority  for  the  project,  because  they  will  be  needed  at  a 
minimum  for  verification  of  the  faster  hybrid  algorithms. 

The  key  elements  of  the  2D  adaptive  grid  RRC  method  [24]  have  been  developed  and 
applied  to  various  problems  [25-26]  in  2D2V  approximation.  To  allow  arbitrary  direction 
of  the  tether  orientation  to  Earth’s  magnetic  field  it  has  to  be  expanded  to  2D3V.  Another 
important  extension  is  an  addition  of  fixed  potential,  current  absorption  and  corrected 
algorithm  for  mesh  refinement  near  the  tethers’  surface.  If  the  potential  will  approach  the 
MV  level,  then  the  relativistic  equations  of  motion  will  have  to  be  added  (at  least  to 
electrons)  to  replace  current  classical  trajectories. 

For  the  VLF  generation  studies  we  are  developing  a  Maxwell  equation  solver.  It  will 
have  to  take  into  account  both  antenna  currents  and  currents  induced  in  the  plasma.  To 
avoid  possible  time  step  limitations  we  will  look  whether  the  so-called  Darwin  model 
[29],  which  in  effect  suppresses  EM-wave  propagation  by  eliminating  transversal 
displacement  currents,  is  beneficial  to  such  studies  from  physical  and  computational 
standpoints. 

Another  avenue,  which  we  are  exploiting,  is  the  introduction  of  an  artificial 
permittivity  of  the  medium  in  the  vast  region  beyond  sheath  and  pre-sheath: 

s  div  Es  =  e^(ft-  fe)dv  (3) 

This  transformation  will  decrease  the  magnitude  of  the  polarization  field,  and  will 
increase  Debye  length  and  thus  reduce  grid  requirements  by  an  order  of  magnitude. 
Unlike  the  sheath,  remote  areas  remain  always  quasi-neutral,  and  charge  separation  does 
not  seem  to  affect  physics  of  the  important  processes  that  occur  in  there.  To  maintain 
undisturbed  Maxwell  plasma  distribution  in  the  periphery  of  a  simulation  domain  we  will 
apply  either  collisonless  technique  [22]  of  continuous  distribution  function  adjustment, 
or,  following  continuous  non-Monte-Carlo  method  [24-26],  we  will  introduce  effective 
collisions. 
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4.1  Features  of  the  Kinetic  Model.  Benchmarking 

The  most  important  pieces  of  model:  particle  trajectory  integrator  and  electro¬ 
magnetic  solver  have  been  developed  in  2D2V  approximation.  In  the  electrostatic  limit 
the  dynamics  of  plasma  species  is  described  with  set  of  kinetic  equations  for  distribution 
function  of  plasma  species  -  electrons  and  ions,  if  neutral  dynamics  is  not  important: 


Of  a  r,  Of  a  |  ( E_EXT  +  ^  +  V  x  B  fxt  )  dfa  _  g 

dt  dx  mac  8w 


a  =  e,i 


(4) 


here  the  self-consistent  electric  field  is  described  by  Poisson’s  equation  depending  on  the 
spatial  distribution  of  charges  as  given  by  Eq.  (3),  while  the  external  electric  (tether 
potential)  and  (Earth’s)  magnetic  fields  are  fixed.  The  right-hand  side  describes  particle 
sources  and  sinks  inside  the  domain  (capture)  and  at  the  boundary  (orbital  motion). 

Trajectory  integration 

For  orbit  integration  instead  of  finite  differences  approximation  of  Newton’s 
equations  of  motion,  we  use  partial  analytical  solutions.  For  example,  in  the  simplest 
2D2V  planar  case  ( B  is  normal  to  velocity)  Eq.  (4)  reduces  to: 


^  +  v  ^  +  v  ^  -(E  +  v  B  )-^- 

8t  x  dx  y  dy  K  x  y  zJdpx 
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immediately  giving  four  equations  for  trajectories: 
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which  could  be  analytically  integrated: 
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B. 


Thus,  we  obtain  an  unconditionally  stable  algorithm,  which  guarantees  that  a  particle  will 
remain  on  a  helical  trajectory  for  arbitrary  time  steps  and  electro-magnetic  field 
magnitudes.  A  few  examples  of  the  calculated  particle  trajectories  in  the  laboratory 
(plasma)  frame  and  in  the  flying  tether  frame  are  collected  in  the  following  Fig.  8. 
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Figure  8  -  The  same  particle  orbit  calculated  in  the  laboratory 
frame  (blue)  and  moving  frame  (red) 

To  correctly  describe  particle  trajectory  in  the  moving  tether  frame  one  has  to  take 
into  account  simultaneously  two  factors: 

i)  particles  have  extra  velocity  V_parallel*  =  V_parallel  +  V_orbit 

ii)  particles  drift  across  B  with  “external”  electric  field  E_normal  =  V_orbit  x  B 

As  follows  from  Fig.  8  these  guarantee  that  the  Larmor  radii  and  velocity  magnitudes 
are  preserved  in  both  systems. 

Another  important  observation  is  the  need  of  actual  electron/ion  mass  ratio  for  the 
predictive  simulations.  In  the  following  two  Fig.9  we  present  calculated  ion  trajectories 
for  commonly  used  artificial  mass  ration  M/m=100  and  for  the  actual  ratio  for  oxygen 
ions. 
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Figure  9  —  Particle  trajectory  and  velocity  plot  (two  upper  panels)  for  artificial 
mass  ratio,  and  (two  lower  panels)  —  the  same  for  actual  oxygen  ions 

As  one  can  see  there  is  a  part  of  a  trajectory  when  ion  moves  in  the  negative  axial 
direction,  breaking  the  strong  meso-thermal  character  of  orbital  motion  (Mi~10). 

Electrostatic  field  solver 

For  Poisson  solver  on  a  structured  mesh  we  deploy  a  fast  FFT-based  solver,  while  an 
SOR-based  solver  is  used  in  the  general  unstructured  grid  case.  In  the  latter  case  the 
benchmarking  problem  was  chosen  to  be  a  field  created  by  a  uniformly  charged  wire.  The 
results  of  comparison  of  the  calculated  and  analytical  potentials  are  given  in  Fig.  10 
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Figure  10-  (left)  -  comparison 
of  calculated  and  analytical 
potential  and  electric  field 
profiles,  and  (right)  — 
exponential  convergence  of 
iterative  procedure 


The  overall  temporal  resolution  and  conservations  properties  of  the  coupled  kinetic 
solver  and  field  solver  are  benchmarked  on  a  plasma  blob  oscillations  in  a  large 
simulation  domain.  Electrons  are  initially  displaced  from  the  heavy  ion  core.  As  follows 
from  Fig.  11  the  period  of  Langmuir  plasma  oscillations  and  the  integral  kinetic  + 
potential  energy  are  maintained  within  1%  accuracy. 
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Figure  11  -  Evolution  of  electrostatic 
energy,  kinetic  energy,  and  spatial 
distributions  of  electrical  charge  and  ion 
density  in  the  test 


4.2  Single  Tether  Simulation 


First  simulations  in  a  large  domain  on  a  structured  grid  involving  the  region  of 
influence  of  a  tether  with  relatively  low  potential  ~~1 00V  indicate  that  i)  the  region  of 
influence  stretches  as  far  as  several  meters  from  the  wire,  ii)  a  separatrix  between  open 
plasma  streamlines  and  trapped  plasma  is  formed,  iii)  there  is  a  large  void  region  formed 
behind  the  wire  and  a  pre-cursor  structure  in  front  of  it.  These  features  are  demonstrated 
by  the  following  Fig.  12. 
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Figure  12  —  Quasi-stationary  flow  around  positively  biased  bare  tether  in  a 
large  simulation  domain 


4.3  Simulation  of  Multiple  Tethers 


First  simulations  of  two  equally  biased  bare  tethers  (on  a  structured  2D  grid  indicate 
the  complex  non-linear  nature  of  the  plasma  flow  interaction.  In  Fig.  13  we  present 
contours  of  electrostatic  potential  (top)  and  ion  density  for  initial  and  steady-state  flows. 
Formation  of  initial  vortices,  density  voids  and  interference  of  the  two  frontal  shock-like 
structures  is  clearly  visible. 


Figure  13  -  (top)  —  contours  of  electrostatic  potential,  and 
(bottom)  -  contours  of  ion  species  density  during  the  initial 
transient  (left)  and  near-steady-state  regimes  (right) 
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4.4  Tethers  with  Negative  Potential 


A  possibility  of  kinetic  modeling  of  a  couple  of  negatively  biased  tethers  that  can  be 
used  for  scattering  of  the  electron  RB  is  presented  in  Fig.  14.  The  formation  of  a  gigantic 
void  in  the  plasma  density  can  be  clearly  seen. 


Figure  14  -  (top)  contour  of 
electrostatic  potential  around 
negatively  biased  tethers  and 
(bottom)  -  of  the  plasma  density 


4.5  Example  of  Adaptive  Grid  Simulation 

Adaptive  grid  simulation  of  a  single  tether  (Fig.  15)  shows  much  finer  detail  of  the 
flow  structure,  which  appears  to  be  more  “turbulent”  compared  to  the  coarser  uniform 
mesh  results  presented  in  Fig.  13.  The  edaptive  branch  is  being  actively  tested. 


I 
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Figure  15  -  contours  of 
plasma  density  obtained  in 
the  adaptive  grid  simulations 
that  automatically  follow 
gradients  of  the  numerical 
solution 
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4.6  Further  developments 

The  focus  of  current  kinetic  model  development  is  the  addition  of  an  arbitrary 
direction  of  the  magnetic  field  with  respect  to  x-y  plane.  This  requires  introduction  of  3D 
velocity  space,  which  substantially  slows  computational  speed,  roughly  by  one  order. 
This  makes  parallelization  of  the  code  using  MPI  [23]  the  next  priority. 

The  following  step  will  include  addition  of  a  Maxwell  set  solver  for  EM-response  to 
the  self-generated  currents  in  plasma: 


1  8E  -  An  -  -  \n  ^  r  _ 

~—  =  rotB - j  =rotB - I  f«vdv 

c  at  c  c  a  J 

- =  —rotE 

c  dt 

divE  =  Anp  =  4/r^q^j  fa  dv 


(8) 


[divB  =  0 

here  index  a  =  e,i  represents  plasma  species.  Neutral  gas  species  can  be  included  into  the 


Ultimately,  3D3V  capability  will  be  very  desirable,  but  will  require  serious  efforts, 
well  beyond  original  plans. 

Extension  to  3D 


However,  two-dimensional  Recursive  Refinement  and  Coarsening  adaptive  grid  was 
expanded  to  3D  dimensions  using  the  following  data  format: 

3D  initial  mesh  (static): 

To  describe  initial  mesh  we  need  two  numbers  -  number  of  basic  elements  = 
elements_ba.se  and  number  of  vertices  =  vertices _base.  Also  we  need  to  specify 
connectivity  list  (6,  elements jbase),  and  provide  coordinates  of  vertices 
x_vertices ( vertices _base) ,  y_vertices ( vertices jbase) ,  z_vertices ( vertices jbase) 

3D  adaptive  mesh  (dynamic): 

elements  -  actual  number  of  elements 

value(V, element)  -  cell-centered  values  of  total  V  unknowns 

level(element)  -  subdivision  level  of  elements,  zero  means  the  element 

belongs  to  the  base  mesh. 

index  (lev el,  element)  -  index  (history)  of  the  element  as  follows 
index(0, element)  -  pointer  to  parent  on  the  base  mesh  Y 

index(0<l<level+l , element)  =  1,. .  .,8  in  accordance  with  3  4 

adopted  sub-cells  numeration:  1  2  X 

neighbor (1 :24,  element)  -  list  of  neighborhood,  with  the  following  7  8 

assumption:  1-4  E,  5-8  W,  9-12  N,  13-16  S,  17-20  Up,  21-24  Down  5  6 
neighbor (1,5, 9,1 3, 1 7,21  ;element)  always  non-zero;  Z 

neighbor(2,6,10,14;18,22;element)=0  means  just  1  in  the  direction 
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neighbor (1,5, 9, 13, 17,21; element)  <  0  means  there  is  a  boundary 
3D  Class- function: 

XYZ(element)  returns  coordinates  of  8  corners  of  the  element  and  of  its  center 

Full  algorithm  includes  local  refinement,  which  was  implemented  and  tested  up  to  3 
levels  of  subdivision,  see  example  in  Fig.  16. 


Figure  16:  Mesh  adapted  to  the  3-D 
current  in  the  thin  wire; 

3  levels  of  subdivision  are  seen  from  2D 
projections  on  X-Y plane; 
also  shown  is  vector  field  of  the 
calculated  B-field 


D 


X 


Current  effort  includes: 

i)  Testing  automatic,  up  to  12  levels,  mesh  refinement; 

ii)  Development  of  the  reversed  function  -  unstructured  mesh  coarsening; 

iii)  Optimization  of  the  data  base  and  two-way  algorithms. 
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These  steps  are  the  cornerstone  of  current  work.  Methodologically  they  are  similar  to  the 
2D  case.  However,  the  complexly  of  the  code  is  one  order  higher.  The  main  reason  is 
algebraic  increase  in  the  number  of  immediate  neighbors,  and  topological  complexity  due 
to  the  need  of  finding  neighbors  of  neighbors  of  neighbors,  one  level  more  than  in  the  2D 
case. 

Another  major  problem  is  continuous  tracking  of  the  particles.  For  efficient  particle 
locating  on  the  unstructured  grid  particles  are  spatially-ordered  with  the  stored  relative 
coordinate  with  respect  to  the  element  they  belong  to.  Tracking  particle  algorithm  is 
substantially  more  complicated  in  three  dimensions.  It  is  being  debugged  in  the  general 
case,  including  the  most  complicated  with  1:8  variations  of  the  ratios  of  linear  elements’ 
sizes  along  cell  comers.  In  2D  case  only  1 :4  corresponding  ratio  was  possible. 


Relativistic  orbits 

The  system  of  Boltzmann  equations  (4)  in  the  relativistic  3D3P  formulation  is  normalized 
using  natural  units  for  the  RBR  problem: 

Elementary  charge:  e  =  1 .6022  x  10  19  C 
Proton  mass:  mp  =  1 .6726  x  10“27  kg 
Length  unit:  1  meter 

Speed  of  light  in  vacuum:  c  -  2.9979  x  108  m  /  sec 

From  these  primary  units  we  derive  secondary  units  for  time,  EM-fields,  kinetic  and 
potential  energies,  etc.  For  instance,  the  dimensionless  equations  of  motion  of  individual 
particles  read: 


Pi 

+  pi 


i  =  \,N 


^  =  qi(E(ri)  +  vixB(ri)) 
l  at 


(9) 


here  pl  =  m0iyivi  is  the  relativistic  momentum  of  the  z'-th  particle,  m0i  is  its  mass  at  rest, 
q(  is  the  integer  (normalized)  charge  and  normalized  relativistic  factor 


Equations  (9)  are  solved  using  semi-analytical  approach  similar  to  [26],  On-going 
benchmarking  includes  calculation  of  orbits  is  analytical  configurations  and  fields 
calculated  by  the  solvers  for  static  EM-fields. 
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Student  involvement 


During  this  project  several  students  received  training. 

Graduate  Student  Christopher  Zeineh  (USA),  M.S.  ’05,  submitted  in  the  end  of  2005 
his  Master  of  Science  thesis  entitled  “ Application  of  an  Electrostatic  High-Voltage 
Tether  to  Radiation  Belt  Remediation 

Undergraduate  Student  Nicholas  Edelman  (USA),  B.S.  ’07  was  involved  through  the 
MIT  UROP  program.  His  report  entitled  “ Research  of  the  Environment  Parameters  for 
Radiation  Belt  Remediation  with  a  High  Voltage  Tether ”  is  attached  to  this  report. 

Visiting  Summer  Student  Jose  Manuel  Zorrilla-Matilla,  (to  graduate  from  Supaero, 
Toulouse  and  from  the  Polytechnical  U.  of  Madrid,  Spain  in  2007)  has  developed  the 
model  for  wave-based  RBR. 

Graduate  student  Yu  Takiguchi,  from  Japan,  has  been  working  on  the  emission  of 
Whistler  waves  by  a  tether,  and  intends  to  submit  this  work  as  his  MS  Thesis  in 
September  2008. 
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APPENDIX  1:  THE  RBR  ENVIRONMENTAL  BACKGROUND 


Nicholas  Edelman,  January  10,  2006 


i)  Environment  Parameters  for  Radiation  Belt 

ii)  Possible  remediation  with  a  High  Voltage  Tether 


1.  Background  Information 

The  earths  magnetic  field  traps  high  energy  charged  particles  in  the  magnetosphere, 
forming  an  inner  proton  belt  and  an  outer  electron  belt.  High  energy  particle  fluxes  lead 
to  radiation  damage  and  degradation  of  spacecraft  and  satellites.  A  high  altitude  nuclear 
detonation  would  likely  create  a  new  high  energy  electron  belt  and  increase  the  radiation 
risk  to  spacecraft  and  satellites.  To  mitigate  the  effects  of  this  new  belt,  a  high  voltage 
tether  system  has  been  suggested  as  a  method  for  accelerating  the  deflection  of  high 
energy  particles  into  their  loss  cone.  On  July  9,  1962,  the  United  States  detonated  1.4 
Megaton  nuclear  warhead  code-named  ’’Starfish”  at  an  altitude  of  400km  (L  =  1.12) 
above  Johnson  Island  in  the  Pacific  Ocean.  The  explosion  released  high  energy  neutrons 
and  fission  fragments  and  ionized  large  quantities  of  neutral  atmospheric  gases.  This 
created  an  artificial  electron  radiation  belt  centered  at  L  =  1.5  [9]  within  the  naturally 
occurring  inner  Van  Allen  proton  belt  (Appendix  II  contains  a  more  detailed  description 
of  the  particles  processes  involved  after  the  explosion).  In  addition,  the  explosion  created 
a  thin  region  of  very  high  energy  electrons  near  L  =  1.2  [6],  On  October  28,  1962,  the 
U.S.S.R.  detonated  a  submegaton  nuclear  weapon  at  L  =  2,  and  as  with  the  Starfish 
explosion,  the  trapped  radiation  belt  formed  close  to  the  L-value  of  detonation  [15],  As 
expected,  the  electron  distribution  best  resembles  the  fission  spectrum  near  the  L-value  of 
detonation  (Figure  1.1).  Although  many  altitudes  are  suitable  for  simulation,  I  propose  a 
study  at  an  altitude  of  3000km  (L  =  1.5)  centered  on  the  equatorial  magnetic  latitude.  The 
lower  atmosphere  (a  few  hundred  kilometers  altitude)  would  not  be  a  suitable  simulation 
zone,  because  the  decay  times  are  in  the  order  of  days  and  thus  are  not  a  primary  concern. 
L  =  1 .5  is  in  a  region  of  substantial  high  energy  proton  flux  and  is  at  a  similar  L-value  to 
the  trapped  radiation  belts  created  by  the  Soviet  and  U.S.  nuclear  detonations. 
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Figure  1.1:  Trapped  electron  spectrum  at  various  L-values  following  the  October 
28, 1962  Soviet  nuclear  detonation  [15]. 


2.  Inner  Van  Allen  Plasma  Environment 

The  inner  Van  Allen  belt  extends  from  approximately  L  =  1.2  to  L  =  2  and  naturally  traps 
high  energy  protons.  In  the  event  of  a  nuclear  detonation,  the  region  would  become 
heavily  populated  with  high  energy  electrons.  For  the  purposes  of  this  analysis,  the  tether 
assumed  to  be  operating  at  L  =  1.5  at  equatorial  magnetic  latitude  (  X=  0),  or 
approximately  an  altitude  of  3000km. 

Table  2.1  summarizes  the  characteristics  of  the  natural  plasma  environment. 
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Characteristics  of  Natural  Plasma  Environment 

Magnetic  Field 

B  =  9.S15xlO-a  Gauss 

Plasms,  temperature  (daytime) 

T  =0.4  eV 

Plasma,  temperature  nighttime 

T  =0.2  eV 

Plasma  density 

n  =  I0*cm-a 

Ion  composition 

H~  =  Sxl0Arm-a 

H?+  =  10acm"a 

0+  =  lOacra-a 

Neutral  density 

n  =  lO^om”3 

Neutral  composition 

H  =  L04cm-a 

H*  =  104cm-a 

O  =  10acm"a 

L armor  Radius 

ra  =  15.3cm 
i[  =  6.57xlt)acm 

Debye  Length 

Ad  =  4.70cm 

Collision ality  (mean  free  path) 

Aa  =  135km 

Aj  =  102km 

Table  2.1:  Environment  Parameters 
2.1  Particle  Energies 

The  dominant  space  plasma  environment  in  inner  Van  Allen  Belt  consists  of  a  relatively 
cool,  low  density  background  plasma  (see  Figures  2.1  and  2.2).  Extrapolating  from  the 
data  in  Figure  2.2,  the  natural  plasma  temperatures  at  an  altitude  of  3000km  and  =0  are 
approximately: 


Plasma  temp-Biaturefdaytime)  =  0.4eV 
Plasma  temper ature( nighttime)  =  0.2eV 


Precipitating  Electrons: 

*  High  Energy:  1 ,0  to  1 00  keV 
■  Low  Density;  1 .0  lo  1 0  om^3 


POUR  -  55°  IO  90°  Latitude 


Increasing 


GEO:  6.6  Earth 
*  High  Energy: 
L  EO:  1 00  to  2 1000  km  *  Low  Density: 

•  Low  Energy:  0.1  to  0.3  eV 

*  High  Density:  1 Q2 lo  1 06  cnr* 


Radii 

1.0  1o  50  keV 
0.1  to  t.O  cm'* 


Figure  2.1:  Properties  of  the  natural  space  plasma  as  reported  in  a  NASA 
report  on  spacecraft  charging  [2]. 
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Figure  2.2:  Variation  of  electron  temperature  with  altitude  and  magnetic  latitude. 

The  data  depicts  the  average  temperatures  from  1989  to  1995  recorded  by  the  thermal 
electron  distribution  instrument  on  the  EXOS-D  satellite  [10]. 


2.1.1  High  Energy  Particle  Fluxes 

In  addition  to  the  average  particle  populations,  the  natural  space  plasma  environment  is 
characterized  by  high  energy  particle  fluxes,  which  are  of  particular  interest  for  studying 
radiation  belt  remediation  applications.  These  fluxes  are  illustrated  in  Figures  2.3  and  2.4. 
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Figure  2.3:  High  energy  proton  fluxes  and  energies  as  a  function  of  L-value 
from  the  AP8MIN  model  [3]. 


2.2  Particle  Densities 

At  the  high  altitude  of  the  inner  Van  Allen  Belt  region,  the  environment  is  characterized 
by  a  relatively  low  plasma  density.  Near  L=1.5  at  an  equatorial  magnetic  latitude,  Figures 
2.5  and  2.6  contain  an  ambient  density  n  =  10  4  cm"3 . 


Figure  2.4:  High  energy  electron  fluxes  and  energies  as  a  function  of  L-value. 
The  data  reflects  Vampola’s  update  of  the  AE-8  model  using  data  from  the 
Combined  Release  and  Radiation  Effects  Satellite  (CRRES)  [3]. 
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Figure  2.5:  Electron  densities  ( m-3 )  described  as  a  function  of  distance  from 
the  earth.  Densities  were  calculated  by  the  Global  Plasmasphere  Ionosphere 
Density  ( GPID )  model  [4]. 
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Figure  2.6:  Variation  of  electron  density  with  altitude  at  the  magnetic  equatorial 
latitude  [12]. 
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2.2.1  Neutral  Densities 


The  neutral  particles  from  Figures  2.7,  2.8(a),  2.8(b),  and  2.9  supports  the  fact  that  the 
high  altitude  space  plasma  is  a  weakly  ionized  plasma  with  neutral  density  approximately 
equal  to  the  plasma  density.  From  the  data,  it  can  be  seen  that  the  neutral  particles  density 
is  nneutrai =  104  cm"3.  The  data  indicates  that  the  neutral  density  composition  is  as  follows: 

H  =  104cm'3 
He  =  104cm'3 
O  =  103cm~3 


2.2.2  Ion  Composition 

From  Figure  2.10,  it  can  be  seen  that  the  inner  Van  Allen  belt  ions  are  dominated  by  low 
energy  protons.  The  region  contains  the  following  ion  densities: 

H+=  5xl03  cm'3 
He  =  103  cm'3 
0+=  102  cm'3 


Figure  2. 7:  Density  profile  of  atomic  hydrogen  with  altitude  as  measured  by 
the  Dynamic  Explorer  ( DE-1 )  satellite  [14]. 
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Figure  2.8:  Density  profile  of  neutral  hydrogen  and  oxygen  during  solar  maximum 
and  solar  minimum  conditions.  The  data  presents  the  output  from  the  Mass- 
Spectrometer-Incoherent-Scatter  (MSIS)  model  [16]. 
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Figure  2.9:  Dependence  of  neutral  helium  and  oxygen  densities  on  L-values 
as  obtained  by  data  from  the  MSIS  model  and  the  IGRF-95  geomagnetic 
field  model  [1 7]. 
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Figure  2.10:  Dependence  of  various  ion  densities  on  altitude  from  the  field-line 
interhemispheric  plasma  (FLIP)  model  [1 7J. 


2.3  Earth’s  Magnetic  Field 


The  earth’s  magnetic  field  can  be  approximated  to  first  order  by  a  magnetic  dipole 
centered  at  the  earth’s  core  and  tilted  1 1°  from  the  rotation  axis. 


B  =  Be 


Re  Vs  [4  —  3  ccs2A] 
~R  )  cob  eA 


(2.1) 


At  R  =  Re+  3000km  and  X=  0,  The  strength  of  earth’s  magnetic  field  is  approximately, 
9.815x10  2  Gauss 
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2.4  Nuclear  Detonation  in  Space  [15] 

Assumption:  A  one  megaton  thermonuclear  bomb  is  detonated  an  altitude  of  a  few 
hundred  kilometers. 

•  Immediately  after  detonation,  massive  radiation  bursts  ionized  nearby  gases.  Because  of 
the  Earth’s  magnetic  field  and  the  particles  in  the  local  plasma,  most  of  this  newly  formed 
plasma  is  trapped  and  reaches  lower  energies  quickly  through  atmospheric  interactions. 

•  A  nuclear  explosion  yields  approximately  1023  neutrons  per  kiloton.  The  neutrons 
ejected  away  from  the  earth  travel  with  such  a  high  velocity  that  most  escape  the 
magnetosphere  before  decaying.  However,  when  neutrons  are  ejected  towards  the  earth, 
many  are  reflected  by  the  atmospheric.  Because  of  interactions  with  atmospheric 
particles,  a  substantial  percentage  of  the  neutrons  decay  in  the  magnetosphere.  Although 
protons  are  produced  during  -  decay,  the  additional  flux  increase  is  negligible  in 
comparison  to  the  ambient  high  energy  proton  flux  of  the  inner  Van  Allen  region. 

•  The  predominant  source  of  particles  arises  from  the  beta  decay  of  fission  fragments.The 
fission  fragments  undergo  -  decay  in  the  following  reaction: 

n  p+  4-  e~  + 

This  reaction  converts  a  neutron  to  a  proton  in  the  nucleus  and  ejects  high  energy 
electrons  and  antineutrinos.  Combined  with  the  decay  of  the  free  neutrons,  this  reaction 
combines  to  create  an  artificial  electron  belt. 

2.4.1  Decay  of  Starfish  Electrons 

After  the  1962  “Starfish”  space  nuclear  detonation,  high  energy  electron  fluxes 
bombarded  the  inner  Van  Allen  region  for  mainly  months.  As  a  result  of  these  high 
energy  electron  fluxes,  many  satellites  in  the  region  suffered  extensive  radiation  damage 
that  led  to  the  failure  of  several  US  and  Russian  satellites.  The  decay  of  the  Starfish 
electrons  is  depicted  in  Figures  2.1 1,  2.12,  and  2.13. 


Figure  2.11:  Average  decay  lifetimes  of  Starfish  electrons  (E  >  0.28MeV ) 
from  the  analysis  of 1963-38C  satellite  data  [9J. 


36 


L  [tarth  radii) 


Figure  2.12:  Decay  lifetimes  (days)  for  Starfish  electrons  ( E  >  0.2MeV )  as 
determined  by  the  data  from  the  1963-3 8C  satellite  [9]. 


Figure  2.13:  Relative  flux  decrease  of  Starfish  electrons  (E  >  5MeV )  with  time  [15] 
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2.5  Sputtering  Yield 


In  the  inner  Van  Allen  radiation  belt,  the  tether  is  subject  to  high  energy  fluxes  of 
particles.  When  high  energy  ions  bombard  the  surface  of  the  material,  the  impact  energy 
causes  the  dislodging  of  surface  particle  from  the  bulk  material.  As  can  be  seen  in  Figures 
3.1(a)  and  3.1(b),  the  sputtering  yield  does  not  increase  monotonically  with  energy.  As 
energy  increases,  the  ions  are  able  to  penetrate  deeper  into  the  bulk  of  the  material  and 
less  of  their  energy  is  transferred  to  atoms  on  the  surface.  Smaller  ions  begin 
experiencing  decay  in  sputtering  yields  at  a  lower  level  than  larger  ions  due  to  their 
smaller  atomic  cross  section. 
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(a)  Dependence  of  sputtering  yield  an  en-  (b)  Sputtering  yield  of  carbon  with  tleu- 
erey  for  aluminum  bombardment  with  hy-  teiium  and  hydrogen  ions.  The  graph 
drogen  ions  [13|.  exhibits  similar  shape  and  magnitude 

to  the  sputtering  of  copper  by  hydrogen 
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Figure  3.1:  Sputtering  yield  for  materials  of  interest  in  electrodynamic  tether 
implementations. 


Materials  suggested  for  use  in  an  electrodynamic  tether  system  include  aluminum  for  its 
high  conductivity  to  mass  ratio  and  carbon  nanotubes  for  their  high  strength.  Figures 
3.1(a)  and  3.1(b)  present  empirical  sputtering  yields  for  aluminum  and  carbon 
bombardment  by  hydrogen.  For  an  upper  bound  calculation  of  the  depth  loss  rate  due  to 
sputtering,  the  following  assumptions  are  made: 


•  Energy  spectrum  consisting  of  IKeV  protons  (the  corresponding  maximum  sputtering 
yield  for  aluminum,  the  weaker  of  the  two  materials) 

•  Proton  flux  of  1010  cm'2  s'1  (higher  than  observed  fluxes  at  any  energy  level) 
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•  For  simplicity,  the  aluminum  is  represented  as  an  infinite  planar  surface. 

Using  the  above  assumptions  the  sputtering  yield  of  aluminum  by  hydrogen  ions  can  be 

Y(E)  &  10-Llbo^fAl 

calculated.  1 .  At  IKeV,  the  sputtering  yield  is  (an  upper  limit 

chosen  from  the  data  in  Figure  3.1(a)).  This  produces  a  loss  rate  of  109  atoms  of  Al/H+, 
or  3.14xl016  atoms  of  A1  per  cm2,  per  year.  Knowing  the  density  of  aluminum  to  be  2.70 
g/cm  ,  the  aluminum  structure  contains  6.026x10  atoms  of  A1  per  cm  .  The  depth  loss 
rate  with  time  in  years  can  be  expressed  as  dL/dt=5.21xlO~5  mm/yr.  Thus,  it  would  take 
1.92x105  years  to  sputter  a  depth  loss  of  1mm. 

Despite  choosing  upper  bound  parameters  to  increase  the  sputtering  yield,  the  sputtering 
of  aluminum  by  hydrogen  did  not  have  a  substantial  impact.  The  ineffectiveness  of  this 
sputtering  clearly  demonstrates  that  sputtering  will  have  a  negligible  impact  on  a  space 
tether  implementation. 
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Abstract 

A  computational  algorithm  is  developed  and  executed  to  calculate  the  rate  of  depletion  of 
magnetospheric  ions  by  an  electrostatic  tether  at  various  altitudes.  This  computation 
relies  upon  past  studies  in  the  OML  regime  of  charged  tethers  to  determine  the  deflection 
angles  incurred  upon  incoming  ions  at  any  given  incident  velocity.  Calculated  depletion 
rates  are  used  to  computationally  estimate  the  time  required  to  deplete  a  given  range  of 
the  magnetosphere. 
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Chapter  1 
Background 


The  Earth’s  magnetosphere  is  defined  as  the  region  within  which  phenomena  are 
dominated  by  the  Earth’s  magnetic  field  and  extends  to  roughly  10  earth  radii  at  the 
equator.  Collisions  among  incoming  cosmic  neutrons  produce  mass  amounts  of  charged 
particles  which  in  turn  have  the  tendency  to  travel  along  magnetic  field  lines,  and 
populate  the  magnetosphere.  Many  such  particles  follow  the  field  lines  into  one  of  the 
Earth’s  poles  and  thus  exit  the  magnetosphere,  but  many  are  repelled  by  the  increasing 
magnetic  field  and  thus  become  trapped  within  the  magnetosphere  for  long  periods  of 
time.  Satellites  passing  through  or  orbiting  within  the  magnetosphere  require  expensive 
shielding  against  high-velocity  trapped  ions,  particularly  those  whose  energies  exist  near 
or  above  the  order  of  1  MeV,  and  we  can  minimize  the  amount  of  damage  incurred  if  we 
artificially  reduce  the  magnetospheric  populations  with  electrostatic  tether  satellites. 
Furthermore,  in  the  event  of  a  nuclear  warhead  detonating  at  high-altitude,  the  saturation 
of  the  magnetosphere  may  be  further  exacerbated,  in  which  case  a  tether  satellite  can 
assist  in  expediting  the  natural  remediation  process. 


1.1  Magnetic  Mirrors 

A  magnetic  mirror  is  defined  as  a  magnetic  field  configuration  in  which  the  field  strength 
changes  along  the  field  lines  in  such  a  way  that  charged  particles  tend  to  reverse  direction 
once  in  the  high-field  region.  Such  a  configuration  usually  consists  of  parallel  magnetic 
field  lines  constricting  and  intensifying  along  a  given  axis.  The  Earth’s  magnetic  field 
produces  a  similar  phenomenon  near  either  magnetic  pole,  near  which  the  Earth’s 
magnetic  field  lines  constrict. 

If  a  charged  particle  is  traveling  along  a  magnetic  field  line  within  a  uniform  magnetic 
field  and  no  electric  field,  then  the  field  exerts  no  forces  upon  the  particle  beyond  those 
contributing  to  Larmor  gyrations.  Thus,  particles  in  uniform  magnetic  fields  free  of 
electric  fields  travel  with  constant  velocity  perpendicular  to  the  field  line  and  with 
constant  radial  speed.  However,  if  the  magnetic  field  lines  begin  to  converge,  a  new 
force  is  introduced  parallel  to  the  field  line. 
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Figure  1-1:  Particle  approaching  a  magnetic  throat. 

Define  a  particle  of  charge  q  and  velocity  v  traveling  along  a  magnetic  field  line  B.  At  a 
given  point  along  the  line,  an  adjacent  field  line  of  similar  magnitude  converges  towards 
it  and  lies  at  distance  r  L.  The  force  on  the  particle  in  the  direction  parallel  to  its  initial 
velocity  is  expressed: 


<  Tf|  >=  — | qx  x  B|  sin  a  =  -|g|v±i?  sin  a 


(1.1) 


sin  a 


In  simplify: 

<  ^1  >=  kh B> 


Calculate  the  gradient  of  B  to  obtain  Br  as  a  function  of  Bz: 

V  B  =-—(rBr)  +  —  B,  =0 

r  or  oz 


Integrate  to  obtain: 


(1.2) 


(1.3) 


(1.4) 
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dB 

Suppose  the  radius  of  the  curvature  of  the  lines  is  small  enough  that  — —  is  constant: 

dz 


M  =-£  rdr 


dB _  _  1  ,  dB 


dz 


2 

2  dz 


Br{rL)  =  -\rLd-^- 
2  OZ 


(1.4) 


Substituting  this  expression  back  into  the  parallel  force  yields: 

1  2 


.  F._  u  v±.rLdBz 


/nv. 


as. 


B  dz 


(1.5) 


We  define  p  as  the  magnetic  moment  such  that 

1  2 

-HIV, 

M  =  (1-6) 

<*ii  >=-/i^  =  //'vnB 

& 


(1.7) 


As  a  particle  enters  a  region  where  the  magnetic  field  lines  converge,  it  experiences  a  net 
parallel  retarding  force.  Depending  on  the  particle’s  velocity  and  the  magnetic  field,  the 
particle  may  be  reflected  back  entirely.  Suppose  that  such  a  particle  does  get  reflected. 
Define  a  series  of  magnetic  field  lines  constricted  on  two  sides.  Now  define  a  particle  at 
the  center,  where  the  field  lines  are  parallel,  as  having  velocity  v0  =  (vin ,  v  0 ) .  At  the 

point  of  deflection,  the  velocity  in  the  direction  parallel  to  the  magnetic  field  is  zero,  so 
v r  =  (v±r,0) .  First  we  employ  energy  conservation: 


\miV\\r2  +VW02)=\mV- 


(1.8) 
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Then  p  conservation: 


1  mv±02  _  1  mvlr2  (1-9) 

2  B()  ~2  Br 


Thus: 


L 

Bn 


=  1  + 


vv<uy 


1 


9()  <  sin  1 


(1.10) 


Thus,  if  the  particle’s  velocity  vector  lies  at  too  great  an  angle  from  the  magnetic  field 
line  along  which  it  travels,  the  particle  will  reverse  direction  and  be  reflected  backwards. 

1.1.1  Magnetic  Mirrors  and  the  Van  Allen  Belts 

Just  as  the  magnetic  field  lines  in  our  mirror  example  expand  and  contract,  so  do  the 
Earth’s  magnetic  field  lines,  which  converge  towards  the  magnetic  poles  and  expand  near 
the  equator.  A  particle  traveling  along  one  of  these  lines  may  have  its  velocity  oriented 
insufficiently  along  the  line  to  overcome  the  mirror  forces  incurred  upon  approaching  a 
pole.  Approximating  symmetry  of  the  earth’s  magnetic  field  lines  across  the  equator,  a 
particle  deflected  close  to  one  pole  will  maintain  its  parallel  speed  and  thus  be  similarly 
deflected  upon  approaching  the  opposite  pole.  Thus,  a  charged  particle  with  insufficient 
speed  to  overcome  the  earth’s  mirror  effect  along  a  given  magnetic  field  line  will  be 
trapped  along  that  field  line. 

The  magnetic  field  lines  along  which  this  phenomenon  is  particularly  prevalent  are  called 
the  Van  Allen  Belts,  or  radiation  belts.  The  high-energy  charged  particles  which  heavily 
populate  the  radiation  belts  prove  hazardous  to  satellites  attempting  to  pass  through  these 
ranges  of  space,  and  sufficiently  protecting  these  satellites  against  the  radiation  is  terribly 
costly,  so  depleting  the  radiation  belts  by  a  given  magnitude  will  cut  costs  on  satellite 
protection  in  future  launches. 

Another  concern  pertaining  to  the  radiation  belts  lies  in  the  possibility  of  nuclear 
intercontinental  ballistic  missiles,  launched  either  in  error  or  intentionally,  being 
intercepted  and  detonated  in  mid-flight.  Nuclear  clean-up  efforts  close  to  the  surface  of 
the  Earth  take  many  years  already,  but  the  radioactive  debris  emitted  to  the  ionosphere 
and  magnetosphere  by  a  high-altitude  nuclear  detonation  is  trapped  in  the  radiation  belts 
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by  the  mirror  effect.  In  such  an  event,  man-made  satellites  could  help  expedite  the 
radiation  belt  remediation,  thus  reducing  any  possible  ecological  ramifications. 

1.1.2  The  Loss  Cone 

For  a  particle  with  a  given  velocity  vector  traveling  along  a  given  magnetic  field  line,  we 
define  the  range  of  velocities  for  which  the  particle  can  escape  as  the  “loss  cone.”  The 
angle  of  this  cone  is  defined  by  the  limits  in  (1.10)  We  illustrate  this  cone  in  velocity 
space.  For  future  reference,  we  establish  axes  for  a  particle  in  the  ionosphere. 

x:  east  (x>0) 
y:  north  (y>0) 

z:  away  from  the  earth,  “up”  (z>0) 

Since  this  loss  cone  applies  to  all  particles  traveling  along  this  magnetic  field  line,  we  can 
assume  that  all  such  particles  whose  velocities  fall  into  the  loss  cone  have  already  exited 
the  radiation  belts.  Our  objective  is  to  deflect  the  remaining  particles  such  that  their  post- 
deflection  velocities  do  fall  into  the  loss  cone,  and  to  that  end  we  employ  a  dual-tether 
satellite. 

1.1.3  Solution  via  Electrostatic  Tethers 

If  we  attach  an  electromagnetic  bare  tether  to  a  satellite  in  the  magnetosphere,  we  can 
induce  a  potential  difference  between  the  tether  and  the  ambient  plasma.  This  would 
induce  an  electrodynamic  force  upon  the  plasma,  including  the  trapped  radiation,  and 
either  collect  the  charged  particles  or  deflect  them  at  various  angles.  Theoretically,  we 
can  use  such  a  tether  to  deflect  trapped  radiation  belt  particles  such  that  their  post- 
deflection  velocities  would  fall  into  the  loss  cone.  Much  research  has  been  conducted  on 
the  behavior  of  tether  satellites  with  regards  to  collecting  charge,  and  much  of  this 
research  can  be  used  to  analyze  its  deflecting  properties  as  well. 

1.2  Previous  Research 

Much  prior  research  has  been  conducted  on  the  behavior  of  electrodynamic  tethers  in 
plasma,  including  computational,  theoretical,  and  experimental.  Since  this  thesis  is 
focused  upon  theoretical  trajectories  of  charged  particles  in  the  vicinity  of  the  tether,  I 
will  review  some  of  the  most  pertinent  theoretical  publications  and  their  most  useful 
equations.  Computational  and  experimental  research  projects  are  also  utilized  to 
establish  the  proper  parameters  for  our  calculations  and  will  be  sourced  as  they  warrant. 

1.2.1  Langmuir  and  Mott-Smith  -  1926 

Langmuir  and  Mott-Smith  pioneered  the  study  of  current  collection  by  spherical  and 
cylindrical  probes,  from  which  they  derived  the  collection  limits  for  thin  and  thick 
cylinders,  named  the  Orbital  Motion  Limit  and  Langmuir  Limit,  respectively.  Electrons 
sense  the  presence  of  the  tether  only  within  an  imaginary  cylindrical  barrier  called  a 
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sheath  (beyond  which  the  potential  difference  is  negligible)  and  their  trajectories  are 
deflected  toward  the  tether  as  they  approach.  In  the  OML  limit,  whether  or  not  it  is 
captured  depends  upon  its  angular  momentum.  In  the  Langmuir  Limit,  the  sheath  can  be 
regarded  as  totally  flat,  so  that  attracted  particles  that  enter  the  sheath  are  all  collected. 
Per  unit  of  collecting  body  area,  the  OML  limit  gives  the  highest  current  possible. 

We  use  this  expression  to  determine  the  tether  architecture. 

1.2.2  Sanmartin  and  Estes  -  1999 


Taking  the  Maxwellian  distribution  for  the  particle  distribution  function,  Sanmartin  and 
Estes  calculated  not  only  the  limit  of  the  probe  radius  for  the  current  collection  to  be 
within  the  OML  regime,  but  also  an  approximation  of  the  collected  current  once  the 
tether  is  within  the  OML  regime.  The  latter  is  used  to  approximate  the  current  collected 
by  the  tether  in  the  magnetosphere,  which  we  use  to  determine  the  tether  architecture. 


I 


OML 


~  2 RLenx 


2e($)  Tether 

V  me 


(1.9) 


1.2.3  B.M.  Minor,  Tethers  Unlimited,  Inc.  -  2003 

Minor  calculated  the  scattering  rates  of  electrons  by  an  electrostatic  tether  composed  of 
several  tethers  bound  together  in  a  cylindrical  alignment,  parallel  to  one  another. 
Employing  analytical  calculations  for  the  two-dimensional  cross-section  of  the  tether 
cluster,  Minor  calculates  the  potential  of  the  tether  and  the  electron  flux  rates  for  various 
energy  levels.  The  electron  results  lie  outside  the  scope  of  this  thesis,  which  will  focus 
solely  on  ion  fluxes,  yet  the  analytical  methods  employed,  specifically  those  pertaining  to 
the  2-dimensional  quantifying  of  tether  potential  and  particle  trajectory,  bear  direct 
application  to  theoretical  electrodynamic  research. 


1.3  Thesis  Outline 

First  and  foremost,  we  will  have  to  determine  what  sort  of  tether  we  wish  to  employ. 
Many  sorts  of  tether  designs  have  been  theorized,  including  single  and  dual  tethers,  yet 
each  of  them  exhibits  its  own  strengths  and  weaknesses  when  faced  with  the  task  of 
deflecting  particles  in  the  magnetosphere.  We  need  to  determine  the  best  design,  or 
perhaps  develop  one  of  our  own. 

Once  we  settle  on  our  tether’s  architecture,  we  will  need  to  determine  the  theoretical 
scattering  properties  of  the  tether  with  respect  to  incoming  particles.  Approximation 
methods  for  extreme  cases  will  also  be  useful. 
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After  determining  the  theoretical  deflection  properties,  we  can  calculate  the  rate  at  which 
the  tether  will  scatter  particles  into  the  loss  cone  when  immersed  in  plasma  with  a  given 
distribution  of  particle  velocities  and  energies.  This  rate  will  determine  the  amount  of 
time  required  for  a  tether  satellite  to  depopulate  a  given  region  of  the  Van  Allen  Belts  by 
a  given  factor. 
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Chapter  2 


Tether  Design 

2.1  Tether  Orientation 

For  a  base  case,  we  assume  an  equatorial,  circular  orbit.  In  order  to  maximize  the  tether’s 
effectiveness,  we  want  to  orient  the  tether  such  that  the  maximum  fraction  of  the  particles 
which  interact  with  the  tether  can  be  redirected  such  that  their  velocities  fall  into  the  loss 
cone.  We  are  utilizing  a  tether  that  is  much  longer  than  its  sheath  is  wide,  so  in  analyzing 
its  electromagnetic  effects  on  a  charged  particle,  we  can  assume  that  the  tether  is 
infinitely  long.  Thus,  a  charged  particle  approaching  the  tether  will  not  undergo  any 
change  in  velocity  in  the  direction  parallel  to  the  tether. 

Next  we  must  determine  the  ideal  orientation  of  the  tether  with  respect  to  the  Earth.  If  we 
orient  the  tether  directly  parallel  to  the  magnetic  field  at  the  equator  (i.e.  in  the  y- 
direction),  the  only  change  in  velocity  will  occur  in  the  vx-vz  plane  of  velocity  space. 
Since  the  tether  is  so  much  more  massive  than  any  incoming  particle,  we  assume  their 
interaction  to  ultimately  result  in  only  a  change  in  the  x-z  component  of  the  particle’s 
momentum,  not  a  change  in  its  magnitude.  Thus,  for  a  particle  interacting  with  a  tether 
oriented  in  the  y-direction,  the  particle’s  velocity  vector  can  only  be  rotated  about  the  y- 
axis  and  not  altered  in  any  other  way.  The  loss  cone  is  radially  symmetric  about  the  y- 
axis,  so  if  a  particle’s  velocity  vector  lies  outside  the  loss  cone  before  interacting  with  the 
y-oriented  tether,  the  velocity  vector  cannot  be  redirected  into  the  loss  cone.  Thus,  a  y- 
oriented  tether  cannot  scatter  trapped  particles  into  the  loss  cone  and  is  inadequate  for  our 
design.  The  tether  must  ideally  lie  anywhere  within  the  xz-plane,  or  equatorial  plane,  for 
maximum  scattering  effect. 

Within  the  equatorial  plane,  the  tether’s  orientation  angle  has  no  effect  on  its  ability  to 
scatter  particles  since  we  are  assuming  the  loss  cone  and  particle  density  to  be  radially 
symmetric.  However,  if  the  tether  is  parallel  to  neither  the  z-axis  nor  the  x-axis,  the 
earth’s  gravitational  field  will  exert  a  greater  force  on  the  end  closer  to  the  earth,  resulting 
in  a  torque  on  the  tether.  Orienting  the  satellite  parallel  to  the  z-axis,  as  proposed  by  TUI, 
thus  results  in  an  unstable  equilibrium.  On  the  other  hand,  if  the  tether  experiences  no 
other  forces  of  similar  or  greater  magnitude,  stable  equilibrium  is  maintained  by  aligning 
the  tether  parallel  to  the  z-axis  so  that  it  points  towards  the  earth,  as  in  Figure  2-1. 
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2.2  Tether  Architecture 
2.2.1  Single-Tether  Designs 

The  simplest  tether  design  is  a  single  bare  tether  attached  to  a  power  supply  such  that  the 
surface  of  the  tether  is  potentially  biased  with  respect  to  its  surrounding  plasma.  Since 
charge  cannot  accumulate  on  the  tether,  a  second  “electrode”  is  needed  to  collect  the 
opposite  polarity  particles,  so  a  “single  tether”  design  is  unphysical. 

If  the  tether  is  biased  positive,  the  excess  charge  can  be  eliminated  by  affixing  a  cathode 
(which  must  itself  be  biased  negatively)  to  disperse  it  out  of  the  end  of  the  tether. 
Unfortunately,  this  produces  a  net  current  within  the  tether,  and  since  we  have  already 
established  that  our  tether  must  lie  perpendicular  to  the  Earth’s  magnetic  field,  the  tether 
succumbs  to  a  Lorentz  force  constantly  accelerating  or  decelerating  the  satellite, 
depending  on  orientation. 
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Figure  2-2:  Single  Bare  Positive  Tether  with  Cathode 

On  the  other  hand,  affixing  a  cathode  either  in  the  center  of  the  tether  or  onto  each  end 
would  split  the  current  into  two  equal  and  opposite  directions.  Rather  than  induce  linear 
acceleration,  the  Lorentz  force  would  now  induce  a  net  torque,  which  could  be  countered 
by  the  gravitational  gradient  between  the  lower  and  higher  cathodes  when  some 
equilibrium  angle  to  the  vertical  is  reached. 


While  such  a  system  would  be  dynamically  stable  up  to  some  current,  attracting  electrons 
with  such  a  high-voltage  tether  produces  a  very  high  current  which  must  be  rejected  by 
the  cathodes.  Sanmartin  and  Estes  give  the  current  collected  by  an  electromagnetic  tether 
as  by  the  following  equations  [1],  where  j  is  an  index  for  each  of  the  N  different  types  of 
ions  present  in  the  radiation  belts. 


Positive  Tether : 


Negative  Tether : 


Ie  =  2  rTen„L 


2eO, 


m„ 
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Ii=Yj2rTen  00^ 
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-2eO, 


mj 


(2.1) 


(2.2) 
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Figure  2-4:  Single  Bare  Positive  Tether,  Central  Cathode  and  Power  Supply 
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Due  to  the  high  voltage  of  the  tether  and  the  low  mass  of  the  electrons,  the  current  from  a 
positive  tether  becomes  strikingly  high.  A  preferable  design  might  utilize  a  negative 
tether,  whose  resulting  current  would  be  reduced  by  two  orders  of  magnitude,  thanks  to 
the  heavier  ions. 

2.2.2  Dual-Tether  Designs 

If  a  simple  tether  is  biased  negative,  positive  ions  would  be  continuously  attracted  to  the 
tether  and  would  remain  affixed  to  the  surface  of  the  tether  until  enough  electrons  are 
attracted  to  neutralize  it.  One  remedy  for  this  is  to  add  a  power  supply  to  bias  one  portion 
of  the  tether  positively  and  sufficiently  so  that  it  can  attract  enough  electrons  to  neutralize 
all  of  the  incoming  ions.  This  effectively  converts  our  single  tether  into  an  anti-parallel 
dual-tether.  However,  the  electrons  collected  by  the  tether  would  enter  from  the  positive 
side,  and  the  collected  ions  would  accumulate  on  the  negative  side.  Thus,  the  collected 
electrons  would  travel  in  a  single  direction  to  neutralize  the  collected  ions,  again  inducing 
a  Lorentz  force  and  producing  in  the  satellite  unwanted  acceleration. 
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Figure  2-5:  Anti-Parallel  Dual  Tether 
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V.V.  Danilov  has  proposed  constructing  a  10-km  parallel  dual-tether  consisting  of  two 
oppositely-charged  tethers  traveling  parallel  to  each  other  and  perpendicular  to  the 
Earth’s  magnetic  field.  His  model  connects  the  two  tethers  at  one  end  by  a  power  supply, 
architecturally  similar  to  the  anti-parallel  tether,  only  bent  in  half. 


10km 


y 


Figure  2-6:  Parallel  Dual  Tether 

As  electrons  reach  the  positive  tether  and  travel  to  the  negative  tether  to  neutralize  its 
attracted  ions,  they  create  a  current  in  each  satellite  of  equal  magnitude  and  opposite 
direction,  producing  a  Lorentz  force  in  each.  These  forces  cause  the  tethers  to  rotate 
about  the  power  supply.  However,  while  this  adjustment  eliminates  the  linear  Lorentz 
acceleration,  it  requires  additional  architecture  to  maintain  functionality.  For  a  tether 
potential  of  1  MeV,  the  sheath  radius  measures  on  the  order  of  102  m,  so  if  the  two  tethers 
are  hooked  onto  the  same  satellite  and  run  parallel  to  each  other,  their  sheaths  would 
intersect,  complicating  the  scattering  model  and  reducing  efficiency.  Keeping  the  two 
tethers  outside  of  each  other's  sheaths  requires  the  addition  of  at  least  one  insulated  beam 
no  shorter  than  the  sum  of  the  two  sheath  radii,  again  on  the  order  of  102  m.  Even  when 
the  tether  sheaths  no  longer  directly  interfere  with  one  another,  their  proximity  dictates 
that  each  tether  would  "block"  a  certain  angular  range  of  incoming  particles  from  the 
other.  This  not  only  restricts  the  angular  range  from  which  each  tether  can  attract  new 
particles  from  the  plasma,  but  it  also  proves  inefficient  as  one  tether  scatters  particles  that 
had  already  been  scattered  by  the  other. 

Furthermore,  as  electrons  reach  the  positive  tether  and  travel  to  the  negative  tether  to 
neutralize  its  attracted  ions,  they  create  a  current  in  each  satellite  of  equal  magnitude  and 
opposite  direction,  producing  a  Lorentz  force  in  each.  These  forces  pull  the  tether  further 
and  further  in  opposite  directions  until  the  entire  tether  satellite  straightens  itself  out  and 
the  Lorentz  force  becomes  linear  one  again.  This  can  be  corrected  by  once  again  adding 
an  insulated  beam  connecting  the  tethers,  this  time  at  the  end  opposite  that  of  the  power 
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supply,  yet  we  continue  to  search  for  a  design  that  is  self-correcting  and  avoids  such 
inefficiencies. 


2.2.3  Series  Tether  Design 


We  thus  propose  a  tether  design  which  we  shall  call  the  “series  tether.”  It  consists  of  a 
series  of  three  tethers  connected  by  two  separate  power  supplies  running  current  in 
opposite  directions.  The  central  10-km  tether  will  be  positively  biased  while  the  two  5- 
km  tethers  on  either  end  will  be  negatively  biased,  the  idea  being  that  each  half  of  the 
positive  tether  will  attract  a  sufficient  number  of  electrons  to  neutralize  the  ions  attracted 
by  the  negative  tether  on  its  own  side.  Each  half  of  the  tether  thus  produces  an  equal 
amount  of  current,  but  in  opposite  directions.  The  induced  Lorentz  forces  can  be 
approximated  as  originating  on  the  center  of  each  half  of  the  total  tether,  and  pointing  in 
opposite  directions,  resulting  in  a  net  torque.  However,  since  two  power  supplies  reside 
10  km  apart  from  each  other,  they  induce  a  substantial  gravitational  gradient  torque  in  the 
direction  opposite  that  of  the  Lorentz  torque.  The  satellite  thus  reaches  a  point  of  stable 
angular  equilibrium  at  which  the  two  torques  negate  each  other. 


Figure  2-7:  Series  Tether 
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2.2.4  Tether  Angle 


For  a  tether  satellite  orbiting  the  earth,  we  define  the  neutral  position  as  when  the  tether 
points  directly  towards  and  away  from  the  earth.  We  define  0  and  <|)  to  be  the  tether’s 
pitch  and  roll,  respectively,  such  that  X=Rcos0  and  Z=Rsin0,  where  R  is  the  radius  from 
the  satellite’s  center  of  mass.  Taking  Rorbit  to  be  the  radius  of  the  satellite’s  orbit 
measured  from  the  center  of  the  Earth,  the  total  gravitational  gradient  torque  vector  is 
given  by: 


Qgg=  3<V 


Q  = 


//  =  GM  Earlh 


3.98xl014 


(2.3) 


Since  our  satellite  will  retain  orbit  within  the  equatorial  plane,  gravity  will  induce  a 
torque  about  only  the  y-axis,  so  the  only  non-zero  term  in  the  gravitational  gradient 
torque  is  the  y-component.  Furthermore,  our  tether  has  no  reason  to  roll,  making  <|)=0, 
and  thus  cos<|)=l .  The  y-component  of  the  torque  is  thus  expressed: 

{QGG)y=3Q.2($(x2  - Z 2 )dmjsin 0 cos 6  (2-4) 


Ignoring  the  weight  of  the  tether  itself  for  now,  we  consider  only  the  torque  resulting 
from  the  power  supplies.  The  integral  in  the  torque  expression  thus  comes  to  represent 
what  we  will  approximate  as  point  masses: 

j(x2  -Z2]dm  =  mx(x2  - Z2)+m2(x2 2  -Z2)  (2-5) 

The  positions  of  the  two  power  supplies  are  represented  by: 

Xx=  R  sin  0  X2  =  -R  sin  6  (2.6) 

Zj  =  R  cos  6  Z2  =  -R  cos  6 


Substituting  these  into  the  mass  integral  expansion  and  assuming  the  masses  of  the  two 
power  supplies  to  be  equal  yields: 
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J  (x2  -  Z 2  )dm  =  m((R  sin  O)2  -  (R  cos  6)2 )+  m((R  sin  6»)2  -  (R  cos  0)2 ) 

(2  7) 

=  2mi?2(sin2  6>-cos2  d)  =  -2mR2  cos 2# 

Substitute  this  back  into  our  torque  equation  (2.4)  to  get: 

{Qgg\  =  3Q.2(-2mR2  cos2<9)sin#cos#  (2.8) 

=  -6mQ2R2  cos 20 sin  20 


Assuming  the  tether  pitch  to  be  sufficiently  small  (0«1),  we  can  approximate  cos20«l 
and  sin20«20.  Furthermore,  the  radius  R  from  the  center  of  the  tether  to  either  power 
supply  is  simply  one  quarter  the  length  of  the  total  tether,  thus  making  our  torque: 


( Qgg)v  ~  12wQ 


f  n\2 


o  =  --mn2i2d 

4 


(2.9) 


Now  we  pair  this  against  the  torque  resulting  form  the  Lorentz  torque.  We  first  focus  on 
either  half  of  the  tether  satellite,  such  that  all  of  the  current  within  our  focus  travels  in  a 
single  direction.  If  we  overestimate  the  half-tether’s  current  to  be  roughly  constant 
throughout  its  length,  the  earth’s  magnetic  field  exerts  the  following  Lorentz  force: 


f  -  UB  sin^ 

F  =  —  lx  B  = - — 

UJ  2 

I  =  current  through  either  power  supply 
i  =  total  length  of  tether 
B  =  magnetic  field 

\| i  =  angle  between  wire  and  magnetic  field 


(2.10) 


Since  we  are  orienting  the  tether  in  the  equatorial  plane,  perpendicular  to  the  magnetic 
field,  v|/=  7i/2,  and  sin \\i  =  1.  Thus,  the  Lorentz  force  on  half  the  tether  is  simplified  to: 

F  RB 
F  = - 

2 

(2.11) 
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The  half-tether  Lorentz  force  can  be  approximated  as  being  enacted  upon  the  half- tether’s 
center,  that  is,  on  its  power  supply.  The  distance  from  the  center  of  the  tether  satellite  to 
either  power  supply  is,  again,  one  quarter  of  the  length  of  the  entire  tether.  Thus,  the 
induced  torque  for  either  half  of  the  tether  is: 


KB  _  I£2B 

- R  — - 

2  8 


(2.12) 


Since  our  satellite  is  both  radially  and  axially  symmetric,  the  torque  on  either  side  is 
equal.  The  sum  of  the  torques  is  thus,  simply: 


U2B 

4 


(2.13) 


At  equilibrium,  the  torques  on  the  satellite  will  cancel  each  other  out,  and  we  are  left  with 
an  expression  for  the  equilibrium  pitch  angle. 

(' Qgg  )y  =  0 

--mQ2£26>  +  R-^-  =  0 

4  4 


3mQ2 

(2.14) 


To  quantify  this  expression,  we  start  by  squaring  our  original  expression  for  Q: 


2  //  _3.98xl014 

/?  3  R  3 

^ orbit  ^ orbit 


(2.15) 


Next,  we  approximate  the  mass  m  of  the  power  supplies  to  be  directly  proportional  to  the 
power  requirements,  accumulating  what  we  will  define  as  a  kilograms  per  watt  of  power 
required.  Approximating  a  to  equal  roughly  20  kilograms  per  kilowatt,  or  0.02 
kilograms  per  watt,  we  deduce: 


m  =  aP  =  aIOtether  =  «/(l06 v)  =  2.0 x  104  •  / 


(2.16) 
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For  altitudes  of  less  than  10  Earth  radii,  the  geomagnetic  field  can  be  approximated  as  a 
dipole  field: 


B  = 


Mo  M  t 

47T  R 


orbit 


(-  2  cos  6,  sin  0, 0) 


(2.17) 


Where  we  adopt  conventional  spherical  coordinates  aligned  with  the  Earth’s  dipole 
moment,  whose  magnitude  is  ME  =8.05x10 22  Am2 .  The  polar  angle  0  at  the  equator 
equals  ji/2,  so  the  magnetic  field  lines  at  the  equator  all  point  straight  north,  as  expected. 
The  magnitude  of  the  magnetic  field  in  the  equatorial  plane  is  thus  expressed  as: 


B  = 


Mo  M_e 

4;r  R 


orbit 


(2.18) 


Substitute  the  previous  expressions  for  Q,  m,  and  B  into  our  angle  equation  to  obtain: 


6  = 


IB 


'  th  Me  ^ 

4*  K,tu  ) 


Mo  M  L 


3mQ2 


3am 


tether 


V  ^  orbit  J 


\27Tjua& 


tether 


6  =  3.38-10  4 /m?  =  0.019° 


(2.19) 


The  equilibrium  angle  of  the  tether  is  not  only  miniscule,  but  also  independent  of  altitude 
and  tether  length.  Thus,  we  can  assert  that  the  Lorentz  torque  on  the  satellite  is  negligible 
and  the  tether  remains  in  stable  equilibrium  parallel  to  the  z-axis,  pointing  almost  directly 
away  from  the  Earth. 

The  only  concern  that  results  from  the  Lorentz  forces  now  is  satellite  deformation 
resulting  from  the  Lorentz  torques  pointing  in  opposite  directions.  This  is  a  structural 
concern  and  lies  outside  the  scope  of  this  thesis,  though  further  study  is  recommended. 
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2.3  Tether  Potential  Magnitude 


According  to  the  tether  current  equations  cited  earlier,  since  the  total  length  of  negative 
tether  equals  the  length  of  its  positive  counterpart,  the  magnitudes  of  the  potentials  in  the 
positive  and  negative  tethers  must  be  unequal,  lest  the  difference  in  mass  between 
electrons  and  ions  produce  a  difference  in  net  current. 

We  assume  that  the  physical  dimensions  (total  length,  radii)  of  both  tethers  are  equal. 
Furthermore,  for  altitudes  above  2000  km,  hydrogen  ions  makes  up  the  vast  majority  of 
the  positive  ions,  so  other  ions  such  as  helium  and  oxygen  can  be  neglected,  simplifying 
the  ion  current  expression  to  a  single  term.  Since  we  wish  to  attract  zero  net  current,  we 
set  the  two  currents  equal  to  each  other  so  they  may  neutralize  one  another. 


K  =  i, 


2rTenxL„  2e°r 


=  2rTen00Lj  - 


m„ 


m, 


Or+=-^cDr-^-^l 
m,  1840 


(2.20) 


(2.21) 


(2.22) 


We  conclude  that  the  potential  magnitude  of  the  negative  tether  is  much  greater  than  that 
of  its  positive  counterpart.  This  makes  sense  because  electrons  travel  much  faster  than 
the  heavier  positive  ions  and  using  similar  tethers  of  equal  and  opposite  potential  would 
result  in  the  capture  of  electrons  at  a  far  greater  rate  than  that  of  ion  capture. 
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2.4  Tether  Radius  and  Sheath  Radius 


2.4.1  Tether  Radius 

Before  we  can  conduct  calculations  using  our  proposed  series  tether,  we  must  confirm 
that  the  tether’s  radius  and  sheath  radius  prove  practical.  The  tether  will  attract  a  current 
proportional  to  the  tether’s  radius,  but  it  will  be  limited  by  the  capabilities  of  the  power 
supplies  which  will  be  driving  the  current.  Let  us  define  a  power  supply  with  an  upper 
limit  of  Pmax  watts.  When  two  are  attached  to  the  satellite,  their  maximum  current 
capacity  can  be  deduced  thus: 


(2.23) 

P  P 

i  _  max  _ _ max 

max_  VT  “10 6  volts 

(2.24) 


Each  power  supply  would  be  responsible  for  transporting  the  current  from  one  half  of  the 
central  positive  tether,  so  that  current  must  be  less  than  Imax.  Applying  our  earlier 
equations  for  the  current  resulting  from  either  half  of  the  positive  tether  length  is  thus 
subjects  the  radius  of  the  tether  to  an  upper  limit.  If  we  assume  the  negative  tether 
potential  to  be  10  MY,  and  if  we  assume  all  positive  particles  are  protons,  then  we 
calculate: 


Ie  L  2e®/ 

—  =  2  rTen-  - 

2  2  \  m„ 


r  2eOr+ 

=  rTenxL, - <  /„ 


m, 


r  < 
rT  ~ 


rT  <90.3^l 


p 

max 

m  P 

e  _  max 

mi  _  ^max 

a 

8 

© 

+ 

2eOr+  enxLQ>T  y 

-2eOr  nxL' ]j 

-2(eOr")3 

(2.25) 


According  to  Figure  2-8,  the  total  number  density  of  the  plasma  within  the  radiation  belts 
numbers  on  the  order  of  1010  m'3  so  if  we  assume  the  presence  of  a  very  strong  power 
supply  capable  of  producing  lOOkW,  the  upper  limit  of  the  tether  radius  becomes 
approximately  1  mm,  a  feasible  radius  if  our  tether  is  made  of  tungsten  steel.  It  should  be 
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noted,  however,  that  according  to  simulations  by  Jean-Marie  Deux  in  Figure  2-9,  the 
current  of  the  orbiting  tether  at  voltages  below  100V  can  add  up  to  twice  the  calculated 
OML  current.  This  differential  appears  to  vanish  as  the  voltage  increases,  yet  stays  fairly 
sizeable  around  the  range  of  550V  covered  by  the  positive  tether.  For  the  sake  of 
conservative  calculations,  we  compensate  for  this  possible  phenomenon  by  halving  the 
tether  radius  to  0.5  mm. 


Log10  Electron  Density  Midday  10  June  1996  160°  West  Magnetic  Longitude 
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Figure  2-8:  Plasma  Density  as  a  Function  of  Altitude  [2] 
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Current  collection  with  orbital  velocity  : 


Figure  2-9:  The  dependence  of  current  collection  on  (positive)  tether  bias  in  the  cases 

of  the  O.M.L.  and  orbital  velocity.  [3] 


2.4.2  Sheath  Radius 

According  to  Choiniere  and  Gilchrist,  the  radius  of  the  tether’s  sheath  is  governed  by  the 
ambient  plasma,  the  potential  bias  on  the  surface  of  the  tether,  and  the  radius  of  the  tether 
[4]: 


2.554 


r  V-325 

r 


J 


(2.26) 


We  already  know  that  our  tether  will  bear  a  potential  bias  of  1  megavolt  on  its  surface,  so 
<Dt  =  106  V.  A.[)  is  the  Debye  length  of  the  plasma  and  depends  upon  the  electron 
temperature  and  density  of  the  plasma,  and  thus  depends  upon  the  altitude  at  which  the 
satellite  will  operate. 


=. 


*o*r. 

e2n„ 


(2.27) 
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The  tether  sheath  thus  depends  upon  electron  temperature,  which  in  turn  depends  upon 
altitude.  The  following  graph  dictates  the  ambient  temperature  of  the  radiation  belts 
below  8000km. 


daytime  nighttime 


Figure  2-10:  Variation  of  electron  temperature  with  altitude  and  magnetic  latitude. 
The  data  depicts  the  average  temperatures  from  1989  to  1995  recorded  by  the  thermal 
electron  distribution  instrument  on  the  EXOS-D  satellite.  [5] 


The  small  temperature  differential  between  the  latitudes  5°  and  15°  indicate  that  we  can 
roughly  approximate  the  temperature  in  the  equatorial  plane  (i.e.  at  0°)  using  the  data 
accumulated  at  5°.  Since  the  temperature  varies  dramatically  between  night  and  day,  we 
approximate  day-average  temperatures  for  several  altitudes,  which  will  in  turn  be  used  to 
approximate  their  corresponding  day-average  sheath  radii.  This  simply  involves 
converting  the  previous  equation  to: 


2.554 


^  r  3 


1.325 


\^D  J 


T  eO)  rr, 

ln^_^  =  0 

rT  kT„ 


(2.28) 


For  a  given  altitude,  we  simply  employ  a  computational  zero-solver  to  find  the  radius 
sheath.  We  need  not  worry  about  the  computer  having  to  choose  between  multiple 
solutions,  since  the  above  equation  increases  monotonically  with  rs. 

Altitude  (km)  Electron  Plasma  Sheath  Radius  Sheath  Radius 
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Temperature 

Debye  Length 

Negative  Tether 

Positive  Tether 

(K) 

(m) 

(m) 

(m) 

2000 

3400 

0.040275 

242.68 

1.2346 

4000 

3900 

0.043134 

234.80 

1.1959 

6000 

5100 

0.049326 

220.11 

1.1237 

8000 

6200 

0.054386 

210.00 

1.0739 

Table  2-1:  Sheath  Radii,  Positive  and  Negative  Tethers 

The  negative  tethers  each  bear  a  sheath  radius  of  order  10  m,  while  the  positive  tether’s 
sheath  is  only  on  the  order  of  lm.  In  either  case,  tether  length  dwarfs  the  sheath  radius  by 
multiple  orders  of  magnitude,  allowing  us  to  approximate  the  geometry  of  the  entire 
sheath  as  a  long  cylinder. 


2.5  Magnetic  Field  Effects 

Our  analysis  of  the  deflection  of  particles  by  the  tether  necessitates  an  analysis  of  both 
electrodynamic  and  magnetic  effects  of  particle  scattering. 

2.5.1  Self-Magnetic  Field 

As  the  tether  collects  current,  the  resulting  magnetic  field  creates  a  series  of  closed  field 
lines  surrounding  the  tether.  The  intensity  of  this  field  reduces  with  distance  until  it 
merges  with  the  Earth’s  locally  open  magnetic  field  lines.  The  planar  projection  of  the 
field  lines  on  the  border  between  the  open  and  closed  field  lines  is  called  the  separatrix. 
Since  the  separatrix  is  not  circular,  Khazanov  et  al  stated  that  the  separatrix  and  the 
circular  induced  field  lines  converge  around  radius  where  they  share  equal  perimeters. 
Thus,  the  condition  for  reduced  current  collection  due  to  self-magnetic  field  effect  is  as 
follows  (in  SI  units)  [6]: 


2/zr.  <  P 


3.23 


separatrix 


cos  a 

rs  =  plasma  sheath  radius 


2ns 0c  B 


“  separatrix 


=  perimeter  of  separatix 


a  =  —  -  (angle  between  tether  and  earth's  magnetic  field) 
B  =  magnitude  of  earth's  magnetic  field 


(2.29) 
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We’ve  already  determined  the  angle  between  the  tether  and  the  earth’s  magnetic  field  to 
be  7t/2,  so  a=0,  and  cosa=l.  The  upper  limit  decreases  with  B,  and  our  sheath  radius 
increases  with  B.  Since  B  increases  with  altitude,  the  inequality  is  most  likely  to  be 
satisfied  at  the  highest  possible  altitude,  which  for  the  scope  of  this  study  is  8000km 
above  the  earth’s  surface.  Substituting  the  magnitude  of  the  magnetic  field  at  this  altitude 
and  the  equation  for  OML  current,  the  inequality  in  (2.29)  can  now  be  quantified: 


f 

2 7irs  <  3.23 

V 


1 

2  7T£0C2B 


\ 

2rTenxL 


<3.23 


rTn„L  -e3  O, 


n1  s0c2B 


2  m: 


<  4.1821  mm 


(2.30) 


The  separatrix  radius  measures  only  about  eight  times  the  radius  of  the  tether  itself, 
which  is  to  be  expected,  given  that  the  total  current  amounts  to  only  a  fraction  of  an 
ampere.  Our  sheath  radius  for  both  the  positive  and  negative  tethers  exceed  this  upper 
limit  by  several  orders  of  magnitude,  so  we  deduce  that  self-magnetic  field  effects  have 
negligible  impact  on  current  reduction  and  can  be  ignored  when  calculating  deflection 
angles. 


2.5.2  Magnetic  Gyrations 


The  Earth’s  magnetic  fields  generate  gyrations  in  the  particles  that  must  be  considered 
when  analyzing  the  electrostatic  effects  of  the  tether,  so  long  as  the  Larmor  radius  (or 
gyroradius)  is  of  a  higher  order  of  magnitude  than  the  radius  of  the  tether  sheath.  The 
Larmor  radius  is  represented  in  CGS  units  thus  [7]: 


r  =  2.38 


K 

B 


cm 


Te  =  electron  temperature  (eV) 
B  =  magnetic  field  (G)) 


(2.31) 


=  1.02-10 


BZ 


cm 
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(2.32) 


T,  =  ion  temperature  (eV) 
B  =  magnetic  field  (G) 

Z  =  charge  state 


m  proton 


Again,  protons  dominate  the  ion  population  above  1500km,  making  them  the  only  ions 
whose  Larmor  radii  we  will  consider,  and  for  which  Z=1  and  y=l .  As  stated  before,  the 
magnetic  field  is  given  by: 


B  = 


Mo  Me 

R 


-Telsa  = 


Mo  M  L 


orbit 


R 


-  ■  1 04  Gauss 


orbit 


(2.33) 


Now  we  can  substitute  (2.33)  into  (2.31)  and  (2.32)  to  calculate  the  values  of  the 
magnetic  field  and  Larmor  radius  for  both  ions  and  elections  at  various  altitudes. 


altitude  (km) 

Te(K) 

Te(eV) 

B(Tesla) 

re(m) 

r;(m) 

2000 

3400 

0.29 

1.37xl0"5 

0.094 

4.024 

4000 

3900 

0.34 

7.22x1  O'6 

0.191 

8.192 

6000 

5100 

0.44 

4.26x1  O'6 

0.371 

15.895 

8000 

6200 

0.53 

2.72x1  O'6 

0.641 

27.467 

Table  2-2:  Average  Temperature  and  Larmor  Radius  in  Ambient  Plasma 

The  thermal  Larmor  radius  of  the  ions  is  considerably  less  than  the  radius  of  the  tether, 
but  since  the  voltage  on  our  tether  is  seven  orders  of  magnitude  larger  than  the  average 
thermal  energy  of  the  particles,  the  particles  greatly  increase  their  energy  by  entering  the 
sheath,  thus  increasing  their  Larmor  radii.  Let  us  assume  that  a  particle  enters  the  sheath 
and  increases  its  thermal  energy  2%  of  the  tether’s  total  potential  bias  (that  is,  20  keV), 
its  Larmor  radius  exceeds  the  radius  of  the  negative  sheath  by  one  order  of  magnitude 
and  the  positive  sheath  by  two  orders. 


Altitude 

B(Tesla) 

re(m) 

0  (m) 

2000 

1.37x1  O'5 

24.52 

1050.94 

4000 

7.22x1  O'6 

46.60 

1997.55 
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6000 

4.26x1  O'6 

79.08 

3389.30 

8000 

2.72x1  O'6 

123.94 

5311.95 

Table  2-3:  Larmor  Radii  at  2%  of  Negative  Tether  Energy 
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Altitude 

B(Tesla) 

2000 

1.37x1  O'5 

4000 

7.22x1  O'6 

6000 

4.26x1  O'6 

8000 

2.72x1  O'6 

re(m) 

r.  (m) 

0.5727 

24.546 

1.0886 

46.655 

1.8471 

79.161 

2.8949 

124.06 

Table  2-4:  Larmor  Radii  at  2%  of  Positive  Tether  Energy 

Thus,  the  effect  of  magnetic  gyrations  in  ions  can  be  neglected  in  determining  the 
deflection  angle  in  both  tethers.  The  same  is  not  true  for  electrons,  however,  so  any 
scattering  calculations  that  neglect  the  effects  of  the  Larmor  radius  will  accurately  apply 
only  to  ions. 
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Chapter  3 
Scattering  Theory 


3.1  Electrostatic  Scattering 

To  calculate  the  total  dispersion  rate  of  the  tether,  we  analyze  the  effects  of  a  single 
electrostatic  tether  on  a  single  incoming  particle.  Defining  the  directional  axes  as  we  did 
when  analyzing  the  loss  cone,  and  assume  a  single  tether  of  uniform  potential  parallel  to 
the  z-axis  such  that  the  origin  lies  in  the  center  of  the  tether.  Next,  assume  a  particle 
barely  inside  the  tether  sheath,  within  the  xy-plane,  whose  trajectory  is  radial  inward. 
Since  the  sheath  of  our  tether  is  considerably  shorter  than  the  total  length  of  any  one 
section  of  the  series  tether,  we  assume  that  the  tether  is  infinitely  long  when  calculating 
the  electrostatic  force  on  the  particle.  Under  this  assumption,  the  tether  is  symmetric  in 
both  the  positive  and  negative  z-directions,  and  thus  exerts  no  force  in  the  z-direction. 
Thus,  we  analyze  the  electrostatic  force  on  the  particle  only  in  the  xy-plane  as  we  would 
a  two-dimensional  Coulomb  collision. 
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Defining  g  as  the  particle’s  speed  in  the  xy-plane  just  as  it  enters  the  particle’s  sheath,  q 
as  its  charge,  r  as  its  radial  distance  from  the  tether,  and  cp  as  angular  position,  the  energy 
and  momentum  equations  defining  a  two-dimensional  Coulombic  interaction  are  as 
follows: 


~^m(r2  +  r2(p2^)+ qtf>(r)  =  ~mg2  S 

mr2(p  =  -mgb 


2  2 

V,  +vy 


(3.1) 


(3.2) 


We  obtain  a  formula  for  the  rate  of  change  in  angle  as  a  function  of  radius  by  rearranging 
the  momentum  equation  as  shown: 


-  mgb  _-gb  (3.3) 

2  2 

mr  r 


Substituting  (3.3)  directly  into  the  energy  equation  (3.1)  produces  an  expression  for  rate 
of  change  in  radius  as  a  function  of  radius. 
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.4-  , 

r  mg 


(3.4) 


Square  both  sides  of  (3.4)  and  divide  by  (3.3)  to  get: 


<P 


dr 

d(p , 
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1  2  2 

r  mg 


(3.5) 


Split  the  differential  terms: 


b  dr 

r2  ~b 2  2qj(r) 

y  r2  mg2 


(3.6) 


72 


To  determine  the  total  change  in  angle  that  results  from  the  Coulomb  collision,  we 
integrate  both  sides.  The  left  side  is  integrated  from  the  initial  angle  to  the  final  angle, 
and  thus  integrates  to  the  total  angle  change.  The  limits  of  the  radius  for  the  integral  on 
the  right  hand  side  are  taken  from  rm  to  Too,  where  rm  is  the  minimum  radial  distance  from 
the  tether  that  the  particle  reaches.  Applying  these  limits  of  integration,  we  get: 


<Pf  , 

A<i>=J</«>=  j-j 

Vi  rm 


dr 


2q<f>(r) 
mg 2 


1- 


2#(fJ  =  0 

mg 2 


(3.7) 


Iisregarding  gravitational  forces,  if  there  were  no  electrostatic  forces  resulting  from  the 
tether,  then  the  particle  would  travel  in  a  straight  line,  and  A(p  would  equal  n.  The  angle 
by  which  the  trajectory  is  deflected  from  the  free  flight  path  is  denoted  %,  and: 


X  =  ?r-2A(p 


(3.8) 


By  the  definition  of  the  tether  sheath,  the  potential  of  the  plasma  outside  of  the  sheath 
radius  rx  is  of  the  order  0.5kTi  «  e®  T ,  and  is  thus  neglected.  We  also  neglect  the 

potential  change  due  to  the  space  charge  of  the  charged  particles  in  transit  through  the 
sheath.  Thus,  we  take  the  potential  to  be  Coulombic  for  r<roo  and  zero  for  r>roo,  and  our 
integral  becomes: 


A  <P=] 


b 


dr 


2q<f){r) 

mg2 


(3.9) 


Now  we  must  calculate  the  potential  function  in  the  Coulombic  region.  We  employ 
LaPlace’s  equation  to  the  electric  potential. 


vV  = 


]_d_ 
r  dr 


d£\ |  i  ay  ,  ay  _  p 


(3.10) 


Assuming  the  potential  field  near  the  tether  to  be  similar  to  that  of  the  field  around  an 
infinite  wire,  the  potential  function  varies  only  with  r,  eliminating  the  second  and  third 
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terms  in  the  differential  equation.  Furthermore,  we  assume  the  tether  to  be  surrounded  by 
macroscopically  neutral  plasma,  so  the  charge  density  in  its  vicinity  is  zero. 


V2  0  =  -— 
r  dr 


( 


d  (  d</> 
r 


dr 


dj)_ 

V  dry 

A 


=  0 


dr 


=  0 


(3.11) 


Thus,  the  content  of  the  derivative  is  constant,  or: 


dd) 
r— ^ 
dr 


=  A 


(j)  =  A  In  r  +  B 


(3.12) 


Our  boundary  conditions  dictate  that  the  potential  at  the  surface  of  the  tether  is  Ot  while 
the  potential  at  the  surface  of  the  sheath  is  0. 


r  =  rT  :  ®T  =  A\nrT  +  B 
r  =  r^  :  0  =  Alnroo+B 


(3.13) 

(3.14) 


Combine  the  two  equations  to  get: 


<f>r  =A  ln-^ 


,  -°r 

A  = - — 


r 

In— 


b=®t  In  To 


r 

In  — 


(f)  =  Ahir  +  B  = 


cl)  _  <T)  In  r 

— -lnr  H - - - ^  =  <fr 


In 


In 


r 

In¬ 


in 


(3.15) 
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(3.16) 


Plug  this  into  our  integral  to  get: 


dr 


2q_ 
mg 1 


r 

In— 


2°r 


r 

In  — 


(3.17) 


The  denominator  of  the  first  integral  contains  many  factors  independent  of  r  which  can  be 
grouped  for  simplicity.  Let  us  define: 


^  <Pr 

mg2  ln?k 

rT 


(3.18) 


Thus,  our  integral  becomes: 


(3.19) 


We  define  77  =  — ,  and  thus  ;;x  =  — ,  and  substitute  into  both  integrals  to  simplify: 
r  r 


I  oo 

A^  =  -J 


dr/ 


drj 


-r]  -Tin  — 

7oO 


i*  li  _  r-.2  _  3  In IL 


•Im 

-V 


dr] 


■  + 


J 


dr/ 


n-  I-77  -Tin  — 

V 


r]_  { VTV 


'Im 

=  \ 


dr]  .  _! 

=  +  sin  77 


I-772  -  Tin  — 
V  ?/x 
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(3.20) 


Substitute  this  into  our  equation  for  the  deflection  angle: 


X  =  7r-Vs.(p  =  ^--2  sin  1  T) 


drj 


-rf  -  Ain  — 

7oO 


(3.21) 


3.2  Approximation  Methods 
3.2.1  Hard-body  approximation 

Fast,  high-energy  ions  are  naturally  capable  of  penetrating  deeper  into  the  sheath  of  the 
positive  tether  than  slower,  lower-energy  ions.  If  an  incoming  ion  possesses  significantly 
less  energy  than  the  tether  possesses  potential,  it  may  be  overwhelmed  by  the  tether’s 
repelling  force  as  soon  almost  immediately  after  entering  the  sheath.  When  such  an  ion  is 
repelled,  it  appears  to  almost  bounce  off  the  edge  of  the  sheath  as  though  it  were  a  hard 
body  collision.  We  can  thus  approximate  our  expression  for  the  deflection  angle  for 
large,  positive  values  of  X  and  expect  to  derive  a  solution  similar  to  that  of  an  elastic  hard 
body  collision  between  a  tiny  object  and  a  cylinder. 

We  start  with  the  definition  for  the  r\m: 


l-//m2  -  Ain—  =  0 


(3.22) 


Now  we  convert  the  logarithm  to  an  exponential  expression  which  we  can  then 
approximate  via  first-order  Taylor  expansion: 


l-7m 


Vm  =7« 


;7« 


1  + 


i-7, 


2  3 


A 


7o=+7o 


1~7» 

A 


(3.23) 


Substituting  this  expression  back  into  the  limit  in  our  integral  yields 
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(3.24) 


X  =  n~ 2 sin"1  7,  -2 | 


'7oo  1  + 


\\-vJ  -Tin  — 


7« 


Now  we  expand  the  term  in  the  radical  of  the  integral’s  denominator  via  a  Taylor 
expansion  around  r|,x: 


Rad  =  l-r/2  -  A\n— =  (Rad)  + 

7o0 

'  -  3 


^  dRad ^ 
dr\ 


(7“7oo ) 


(3.25) 


=  1-7.  + 


-27oo  (7-7oo) 

7.  J 


Assuming  A,»l,  we  can  approximate: 


A 

7oo 


(3.26) 


This  gives  us  the  following  approximation  for  the  radical: 


Rad  =  1  -  TjJ 


-  —  (7-7oo) 

700 


_l  p7oo  1+ - ' 

Z  =  ^  -  2  sin  7M-2  l  * 

J<7» 


drj 


,  l-7oo2  -  —  (7-7oo) 

V  7oo 


(3.27) 


To  solve,  we  substitute: 
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R  =  l-f]J  -—{17-TjJ  dR  =  —dr] 

V*  T/oo 


V  =  V, 


f  !  2  A 

1+1-" 


V 


1 


J 


X  =  n-  2 sin-1  r/x  -  2  J° 

=  ;r-2sin_1  77  -2—  f 

'*  T  Jo  V* 

=  ^-2sin-1  4^/l 


->R  =  0  r/  =  rix  ->R  =  l-tjK 
V^dR_ 

7x  r1-^2  <*R 


/i 


(3.28) 


Since  we  are  still  assuming  that  A>>1,  then  >3I«1 ,  and 


h  = 


/ 

roo  COS 

V 


(3.29) 


This  is  the  equation  for  a  particle  colliding  with  a  hard  cylinder,  just  as  we  expected. 


3.2.2  Soft-body  approximation 

On  the  other  side  of  the  coin,  if  an  ion  possesses  a  very  large  amount  energy  in 
comparison  to  the  tether’s  potential,  such  that  A,«l,  it  can  pass  almost  straight  through 
the  sheath  with  only  minimal  influence  from  the  tether.  Such  an  approximation  applies  to 
only  a  small  minority  of  the  total  number  of  ions,  but  it  can  apply  to  both  attracting  and 
repelling  tethers. 
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Since  the  effect  of  the  tether  will  be  miniscule,  we  can  approximate  the  entrance  and  exit 
angles  of  the  ion  to  be  nearly  the  same. 


8  MAX  =  cos~‘  — 
r 


(3.30) 


We  use  this  to  define  the  range  of  the  particle’s  radius  from  the  tether  and  its  position 
along  the  x-axis  with  respect  to  b  and  0: 


r  = 


cos  0 
x  =  r  sin  6  =  b  tan  0 

b 


0 MAX  <6<  0 MAX 


dx  = 


cos  0 


-dO 


(3.31) 


(3.32) 


(3.33) 

We  define  g  as  the  ion’s  velocity  component  in  the  plane  of  the  tether  cross-section.  If 
we  assume  that  the  electrostatic  force  is  so  weak  that  it  accelerates  the  particle  mostly  in 
the  direction  perpendicular  to  the  particle’s  initial  velocity,  we  can  approximate  g  to  be 
constant,  which  we  use  to  approximating  change  in  time  as  directly  proportional  to 
change  in  the  x-direction.  Further  substitution  yields  change  in  time  with  angle. 
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(3.34) 


b  dd 
g  cos2  6 


The  electric  field  within  the  sheath  is  given  by: 


E  =  - 


dcf) 

dr 


d_ 

dr 


{  \ 

r 

In  — 


r 

In  — 
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d>r  cos# 


r  r  v 

In  — r  In  — 


'  T  J 


(3.35) 


Next  we  define  the  force  on  the  particle  in  the  normal  direction  to  obtain  a  differential 
equation  for  the  ion’s  normal  velocity: 


Fn  =  maN  =  in¬ 


dy  „ 


dt 


=  eE  cos  0 


(3.36) 


Substituting  our  prior  equations  produces  a 


,  eE  ..  cos 2  6 

dvN  =  — cos  Qdt  = - — 


m 


m\n 


b  geos2  0 


dO 


dVN 

g 


(3.37) 


Integrating  over  the  total  trajectory  of  the  particle  within  the  sheath,  we  deduce  the  total 
change  in  the  normal  velocity  (that  is,  the  final  normal  velocity)  in  terms  of  the  total 
change  in  angle. 


rA|,w  dvN 

Jo 


g 


^L  =  A0max=A  cos"1 


80 


(3.38) 


Since  the  ion  originally  has  no  normal  velocity,  the  deflection  angle  is  simply  the  angle  of 
the  final  velocity  vector,  or: 


AvV 

x  =  — 

g 


A  cos-1 


(3.39) 


This  approximation  is  only  valid  for  very  low  values  of  X,  and  since  the  result  itself  must 
be  not  much  greater  than  X. 
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Chapter  4 

Numerical  Implementation 


4.1  Flux  Integral 

Since  we  are  primarily  interested  in  the  scattering  of  high-energy  particles,  we  shall 
restrict  our  focus  to  the  high-energy  populations,  whose  mean  energy  is  1  MeV.  Low- 
energy  particles  are  collected  and  scattered  by  the  tether  as  well,  but  their  primary  effect 
lies  in  their  collection  rates  which  in  turn  determine  the  tether  and  sheath  radii.  Since 
such  effects  have  already  been  calculated  in  the  previous  chapters,  the  low-energy 
particles  can  now  be  regarded  as  a  separate  population  and  ignored  as  we  calculate  the 
scattering  rates  of  the  high-energy  particles  at  mean  temperature  T  =  106T 1600  K 

4.1.1  Distribution  Function 

Now  that  we  have  a  computationally  feasible  expression  for  the  deflection  angle  of  a 
given  particle,  we  will  implement  it  to  determine  what  percentage  of  incoming  particles 
will  be  depleted.  A  tether  in  orbit  will  be  attracting  ions  from  many  different  directions 
and  velocities,  and  using  the  equation  for  the  deflection  angle,  we  can  determine  the 
percentage  of  these  ions  that  are  properly  deflected  into  the  loss  cone.  However,  before 
we  can  do  anything,  we  must  first  determine  what  distribution  the  ion  velocities  will 
obey. 

4.1. 1.1  Loss  Cone  Correction 

We  start  with  a  Maxwellian  distribution  of  particle  velocities  in  the  magnetosphere 
plasma,  dependent  only  upon  the  velocity  magnitude  and  normalized  to  integrate  over  all 
velocities  to  the  density  n. 


W 


e 


2nkT) 


m(v2  +v2  +vz2) 
2kT 


(4.1) 


in*',,  z)dvxdvydvz 

All  Velocities 


/•V  =00  pv  =00  pv  =00  /  \ 

j  {  {  /(v* ,  V  ,  vz  )dvx dv  dvz  = 

Jv  =— 00  Jv, ,=— 00  Jv  =— 00  ^ 


(4.2) 


However,  we  are  assuming  that  all  charged  particles  whose  velocities  fall  into  the  loss 
cone  exit  the  magnetosphere,  so  our  distribution  function  must  not  include  such  ions. 
This  exclusion  depends  solely  upon  the  direction  of  the  ions’  velocity  vectors,  not  their 
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magnitude,  while  the  Maxwellian  distribution  function  depends  upon  velocity  magnitude, 
not  direction.  Thus,  if  we  exclude  a  certain  fraction  of  velocity  directions  from  our 
calculation,  then  the  contribution  of  each  velocity  magnitude  to  the  normalizing  integral 
will  be  reduced  by  that  same  fraction,  as  will  the  integral’s  total  value.  To  determine  the 
fraction  by  which  the  normalizing  integral  is  reduced,  we  simply  determine  the  solid 
angle  fraction  not  encompassed  by  the  loss  cone.  Remembering  that  the  loss  cone  is 
projected  in  both  directions  along  the  y-axis,  we  calculate  this  fraction  k  of  the  total  solid 
angle  47i  to  be: 


Y  =  l-£W  =  i_2'24l-cosoJ  =  cosg 
Ak  An 


(4.3) 


If  we  restrict  the  limits  of  our  normalizing  integral  to  only  those  velocity  vectors  whose 
directions  lie  outside  the  loss  cone,  yet  maintain  our  original  Maxwellian  distribution 
function,  our  result  is  reduced  by  a  factor  of  k: 


JJJ  fk  >  vy  >  vz  ]dvxdvydv_  =  n  k  =  n  cos  aL 

Outside 
Loss  Cone 


(4.4) 


Since  we  want  the  normalizing  integral  to  equal  n,  we  divide  the  distribution  function  by 
the  extra  factor,  thus  converting  the  normal  Maxwellian  distribution  into  a  normal  loss- 
cone  distribution: 


cos  aL  \  lnkT 


3  (  2  ,  2  ,  2s 

\-  m(yx  +v  +vz  ) 

m  \  2 - - - 

m  1  e  2 kT 


(4.5) 


JJJ  / lc  (l  >  V  ^ )dvxdvydvz  =  n 

Outside 
Loss  Cone 


(4.6) 


This  makes  sense  because  constricting  the  limits  of  the  normalization  integral  reduces  the 
result,  and  we  compensate  by  increasing  the  distribution  function. 
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4.1. 1.2  Change  to  Cylindrical  Coordinates 

Our  two-dimensional  calculations  on  scattering  theory  were  derived  in  polar  coordinates 
rather  than  Cartesian  coordinates,  so  it  makes  sense  to  convert  our  three-dimensional 
distribution  function  from  Cartesian  coordinates  to  cylindrical  coordinates. 

We  split  up  the  velocity  vector  into  three  components:  the  velocity  component  parallel  to 
the  tether  (vz),  the  velocity  component  within  the  perpendicular  xy-plane  (g),  and  the 
angle  at  which  the  latter  component  lies  from  the  x-axis  (0).  In  other  words: 

vx=gcos0  g  =  Vv/  +  v>’2 

v  =  g  sin  6  0  =  tan-1  — 

y  v 

x  (4.7) 

We  similarly  convert  the  differential  terms  in  the  integral: 

dvxdvy  =  gdgdO  (4.8) 


Our  loss-cone  distribution  is  now  defined: 


/lc(v  z ’§>&)  = 


cos aL  \27ikT 


m 


m(g  +vz  ) 
2  kT 


Jjj §  ■  fLc (vz’g’ e)d8d edv,  =  n 

Outside 
Loss  Cone 


(4.9) 


(4.10) 


4.1. 1.3  Flux  and  Impact  Parameter 

Now  that  we  have  our  distribution  function,  we  can  calculate  the  particle  flux  per  square 
meter  per  second  through  a  surface  area  in  the  magnetosphere’s  plasma,  assuming  the 
radial  direction  to  be  perpendicular  to  the  surface. 

r=  \\\g2  ■  fLAv^g^)dgdedv zm~2s~X 

Outside 

Loss  Cone  id  in 
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To  calculate  the  actual  number  of  particles  passing  through,  we  have  to  determine  the 
area  through  which  the  ions  are  passing.  However,  we  cannot  simply  multiply  the  per- 
square-meter  flux  equation  by  the  area  of  the  tether  sheath  because  not  all  particles  which 
penetrate  the  sheath  do  so  at  a  perfect  right  angle.  This  problem  can  be  circumvented  by 
working  in  terms  of  the  impact  parameter  b. 

Starting  with  a  cross-sectional  diagram  of  the  tether  sheath,  define  a  velocity  vector  v 
such  that  g  is  greater  than  zero  and  0  any  single  value.  Isolate  all  incoming  particles  with 
said  velocity.  Define  a  Gaussian  surface  directly  between  these  incoming  particles  and 
the  tether  sheath,  as  shown,  such  that  the  ram  end  is  flat,  precisely  the  length  of  the  sheath 
diameter,  and  perpendicular  to  the  radial  direction. 


Since  there  are  no  sources  or  sinks  of  ions,  the  flux  of  ions  of  this  given  velocity  out  of 
the  surface’s  rounded  end  (and  into  the  sheath)  must  equal  the  flux  into  the  plate  from 
outside,  which  we  can  much  more  easily  calculate.  The  incoming  particles  have  no 
impetus  to  pass  through  any  portion  of  the  plate  any  more  than  any  other  part  because  the 
flat-plate  portion  of  the  surface  exists  outside  the  tether  sheath,  resulting  in  equal 
distribution  across  this  plate.  For  ions  of  any  given  velocity  vector  with  a  non-zero  radial 
component,  the  flux  area  is  simply  twice  the  sheath  radius  multiplied  by  the  length  of  the 
tether.  This  holds  true  for  all  such  velocity  vectors,  so  the  entire  flux  integral  is  also 
multiplied  by  this  factor,  as  shown: 

rM,„=  JJJ  g,-fu:(v-_,g,e)-2RJ^-dgdm.s-1  <4-12) 

Outside 
Loss  Cone 


However,  when  we  calculate  the  rate  at  which  particles  are  scattered  into  the  loss  cone, 
not  all  impact  parameters  result  in  the  incident  ion  hitting  its  target.  This  is  not  only  a 
limiting  factor  to  the  range  of  impact  parameters  which  we  sum  into  the  integral  but  also 
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a  function  of  the  velocity  vector  itself.  This  term  for  the  impact  parameter’s  acceptable 
width,  which  we  shall  call  Wb,  must  be  completely  nested  within  the  flux  integral. 

rscatter=  JJJ^ '  '  f LC  (Vz  ,g^\ '  Wb  (Vz  '  LT  '  dgdBdv  zS~l  (4-13) 

Outside 
Loss  Cone 


4.1.2  Scattering  Limits 

Before  we  leave  the  integral  as  is,  we  must  realize  that  Wb  is  not  non-zero  for  all  velocity 
vectors.  That  is  particles  incident  at  certain  velocities  have  no  chance  at  being  scattered 
into  the  loss  cone.  To  minimize  the  computational  requirements  for  this  calculation, 
there  is  no  sense  in  spending  processing  time  towards  ions  which  do  not  contribute 
toward  the  scattering  flux.  We  continue  with  our  two-dimensional  plate  diagram, 
identifying  and  excluding  all  such  “hopeless”  ions  from  the  integral  limits  to  manageable 
levels.  Our  goal  is  to  isolate  the  permutations  of  initial  velocities  and  impact  parameters 
such  that  any  particle  bearing  those  initial  conditions  upon  entering  the  tether  sheath  will 
be  scattered  into  the  loss  cone  upon  exiting  the  sheath. 

4.1.2. 1  Axial  Velocity  Component  vz 

We  start  with  the  vz  parameter  and  allow  it  to  take  any  value.  Placing  no  restrictions 
leaves  the  limits  of  vz  the  same  as  in  the  normalizing  integral: 

—  00  <  V,  <00  (4.14) 


To  shorten  our  calculations,  we  remember  that  we  are  assuming  the  particles  to  act  as 
though  our  tether  were  infinitely  long,  thus  resulting  in  the  same  deflection  in  the 
equatorial  plane  regardless  of  the  sign  of  the  z-velocity  component.  Thus,  we  cut  the 
limits  to  our  integral  in  half  due  to  symmetry: 

0  <  vz  <  co  (416) 


This  change  is  accommodated  by  multiplying  the  total  flux  integral  by  2. 

4.1.2.2  Radial  Velocity  Component  g 

Next  up  is  the  g  parameter.  Remember  that  for  a  particle  with  a  given  velocity,  the 
scattering  tether  can  alter  the  particle’s  velocity  vector  only  by  rotating  it  about  the 
tether’s  parallel  axis.  Thus,  if  vz  is  sufficiently  large  compared  to  g,  then  the  total  vector 
cannot  be  rotated  into  the  loss  cone  no  matter  how  it  is  scattered: 
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Thus,  for  a  given  value  of  vz,  g  much  be  at  least  sufficiently  large  that,  if  the  velocity 
vector  were  to  exist  entirely  within  the  y-z  plane,  it  would  lie  parallel  to  the  edge  of  the 
cone.  Larger  values  of  g  place  the  vector  further  inside  the  loss  cone,  while  smaller 
values  of  g  cause  it  to  exit  the  loss  cone.  Thus,  only  values  of  g  within  this  limit  have  an 
opportunity  to  enter  the  loss  cone,  and  no  other  values  of  g  do  not  factor  into  our 
scattering  flux.  That  is: 


- £ —  <  g  <  oo 

tano^ 


(4.17) 


4.1.2.3  Incident  angle  0 

Next,  we  analyze  the  incoming  particle  angle  0  with  respect  to  both  vz  and  g.  If  we 
assume  that  the  ions  will  be  approaching  the  tether  sheath  independent  of  incident  angle, 
then  the  behavior  of  the  particles  is  symmetric  across  both  the  x-axis  and  y-axis.  To 
simplify  our  calculations,  we  multiply  the  flux  integral  by  four  and  limit  our  integral  to 
one-quarter  of  the  total  range  of  0: 


O<0  <  — 

2  (4.19) 


These  limits  must  be  further  constricted  so  as  not  to  include  any  incoming  particles 
whose  velocities  fall  into  the  loss  cone.  The  condition  for  an  ion’s  velocity  vector  to  lie 
outside  a  loss  cone  is: 


tan  aL 


< 


(4.20) 


Since  we  are  assuming  that  no  incoming  particles  lie  within  the  loss  cone,  we  substitute 
the  initial  polar  velocity  components  into  the  above  equation  and  deduce  which  values  of 
0  satisfy  and  thus  may  be  excluded  from  our  calculations. 


tan2  aL  <  -  +2Vz 


v 


y 


g 2  cos2  0  +  vz 
g2  sin 2  0 


g2  sin2  #tan2  aL  <  g2  cos2  0  +  v2  =  v2  -  g2  sin2  0 
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1 


g2  sin2  o{[  +  tan2  aL)=  g2  sin2  6 


vc°s  ar  J 


<  v1 


=±sin" 


vcosa. 


v  8  j 


(4.21) 


Since  we  restrict  0  to  the  first  quadrant,  the  integral  limit  becomes: 


0  <  6  <  sin 


^vcos aL  ' 
8  , 


(4.22) 
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4.1.2.4  Impact  Parameter  b 

Now,  given  an  ion  with  velocity  vector  components  vz,  g,  0,  we  must  determine  which 
values  of  b  allow  for  loss  scattering.  This  is  accomplished  by  determining  which  values 
of  i  produce  such  results,  from  which  we  determine  the  corresponding  values  of  b.  As 
displayed  in  the  following  figure,  we  define  %  such  that  a  positive  deflection  angle 
corresponds  to  a  deflected  particle  whether  the  impact  parameter  is  positive  or  negative: 


Figure  4-2:  Deflection  Angles  for  Positive  and  Negative  Impact  Parameter 

We  first  identify  the  initial  velocity  components  in  terms  of  vz,  g,  0. 

vx  =  -g  cos  & 
vy  =  -g  sin  3 
Vz=v2 

(4.23) 


We  then  define  the  exit  velocity  components  in  terms  of  the  initial  velocity  and  the 
deflection  angle: 


vx  =gcos(6»  +  ^--^)  =  -gcos(6»- x) 

f 

vy  =gsin(6>  +  ;r- j)  =  gsin(<9- %) 

t 


(4.23) 
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V 


*  =  gcos(0  +  7T  + x)  =  -gcos(0  + z) 

f 

vy  =gsin(0  +  x  +  j)  =  gsin(6»  +  %) 

v2'  =  vz  (4.23) 


In  other  words: 

r 

vx  =-gcos(0  +  z) 

=gsin(6»  +  j) 

v/  =  vz  (4.23) 

_  /-for  b  >  0 
\+  for  b  <  0 


Next,  we  plug  the  components  of  the  exit  velocity  vector  into  the  loss-cone  condition. 


tan  aL  > 


4 


v  _,2+v  ’2 


> 


g2  cos2(6>  + j)  +  vz2 
g2  sin2  (0  +  z) 


g2  sin2(#+  j)tan2  aL  >  g2  cos2(#  +  %)+  vz2 

g2sin2(^  +  ^)(l  +  tan2«I)>g2  +vz2 

sin2(6>  +  z)>x  |  vz2 
cos2  aL  g2 


sin  '-(e+x)>  cos2  aL 


1  +  - 


8 


(4.24) 


Isolate  the  angular  terms  to  obtain  our  limit  for  %  in  terms  of  0: 
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\+  for  b  <  0 
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(4.25) 


From  here  we  can  determine  which  values  of  %  result  in  an  ion  of  velocity  vector  v  being 
scattered.  After  that,  we  can  computationally  determine  which  values  of  b  correspond  to 
each  deflection  angle  limit,  and  the  sum  Wb  would  equal  the  difference  between  each  pair 
of  impact  parameter  limits. 


Figure  4-3:  Deflection  Angles  Dependant  upon  Impact  Parameter 

That  is,  if  each  of  the  limit  deflection  angles  %\,  %2,  %3,  %4,  and  determining  their 
corresponding  impact  parameters  bi,  b2,  b3,  b^ 

Wb  (Vz  ,g,0)=  \bl  ~  b2  I  +  1*3  -  b4 1  ^4’26) 
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4.1.3 


Final  Flux  Equation 


Going  back  to  equation  3.1,  we  establish  the  distribution  function  by  multiplying  the 
interior  term.  After  applying  all  of  the  above  limits  to  our  flux  equation  this  turns  into  the 
following  for  our  scattering  flux: 


r 


scatter 


vcosct^ 

g  V 


wb  (vz  ,g,0)-LT  ■  d6dgdvzs  1 


(4.26) 


4.2  Computational  Approximation 


We  have  already  derived  an  expression  which,  for  a  given  velocity,  accepts  an  impact 
parameter  and  returns  the  deflection  angle.  For  this  integral,  however,  our  task  is  to 
accept  a  deflection  angle  and  return  an  impact  parameter.  The  deflection  angle 
expression  cannot  be  simplified  analytically,  so  instead  we  will  develop  a  reference  table 
from  which  we  can  interpolate  values. 

We  construct  a  2-dimensional  matrix  such  that  each  column  index  corresponds  to  a  given 
value  of  b,  and  each  row  index  corresponds  to  a  given  value  of  X.  The  matrix  elements 
themselves  are  the  values  of  x  corresponding  to  the  values  of  b  and  X  corresponding  to 
that  element’s  indices.  In  the  deepest  layer  of  the  flux  integral,  we  can  determine  the 
limits  of  x  as  defined  by  vz,  g,  and  0,  and  use  the  reference  matrix  to  determine  the 
corresponding  values  of  b. 

The  integral  in  this  expression  cannot  be  determined  analytically,  so  for  the  purpose  of 
our  calculations,  we  will  be  employing  quadratic  approximations  via  Matlab.  One 
problem  that  arises  is  that  quadratic  integral  approximations  substitute  the  integration 
limits  directly  into  the  expressions;  doing  so  with  the  upper  limit  produces  an  infinite 
number  since,  by  our  definition  of  rm: 

_ 1 _ _ _ 1 _  1 

-/lln—  Jl-— ’  -,lln—  ^  (4.27) 

V  ^  \  rm 


This  singularity  at  r\m  produces  an  error  which  we  resolve  by  splitting  the  integral  into 
two  parts:  the  main  body,  and  an  addendum  to  approximate  the  area  closest  to  the 
singularity. 

Define  s  as  a  value  of  r\  slightly  smaller  than  r\m.  We  split  up  the  integral  thus: 


92 


(4.28) 


The  first  expression  on  the  right-hand  side  can  be  derived  computationally  because 
neither  limit  encompasses  the  singularity.  The  second  expression  can  be  approximated 
analytically  by  introducing  a  change  of  variables. 

rf='nm-Ti,  0  <  rf<  s,  s « rjm  (4.29) 

From  here  we  expand  the  term  within  the  denominator’s  square  root. 

l-rj2  -Tin—  =  \-{rjm -Tin 7w~77 

>7,  7oo 

(  » V 

=  —  1-  — 

y  Vm)_ 

r  \  c  ,  x 

=  —  -Ain  1--5- 

\^I co  J  V  ) 

f  \  f  ,  \ 

=  l~7m2  -  Tin  ^  +2^’-7'2-Tln  1-^- 

^ _  V^ooy  V  y 

=0 

=  ^mn'-Ti'1  -Tin  1 - 

V  J 


(4.30) 


We  approximate  this  expression  with  a  second-order  logarithmic  Taylor  series. 
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f  A  ^ 

+  ?7' 

l27m  J 

Urn  J 

(4.31) 


Next  we  substitute  this  approximate  expression  into  the  addendum  integral. 

'r  ell]  nr  dr/ 


2  01  „  7m  nm-s 


n"~£  X-r/1  -Ain  — 

V  7« 


A  n 

f  A  ^ 

— ^-1 

+  77' 

—  +  27m 

1 

y 

l27»  J 

A  J 

(4.32) 


From  our  earlier  definition,  r/'=r/m-r/ ,  we  determine  that  dr/'=-dr/.  Furthermore, 

while  the  limits  of  r|  are  r|m-£  and  r|m  ,  the  corresponding  limits  of  rf  are  s  and  0, 
respectively.  Substituting  all  of  this  into  our  integral  yields: 
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■dr/' 
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A  , 

fA  „  1 

77'2 
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+  77' 
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0  2  1 
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1  +  77 


,  1  ^-2  7„ 
27m  /i  +  277n 


(4.33) 


To  make  the  integral  a  bit  more  feasible,  we  compute  the  first-order  Taylor  expansion  of 
the  term  about  77'  =  0 . 


I  +  77' 


1  *-2 A  A 

27m  A  +  277m2y 


«1- 


47m  A  +  27m2 


r/' 


(4.34) 


Plug  this  into  the  integral  to  get 
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dr]} 


(4.35) 


A  ^  o  I — 7  I  1  A  —  2r/ 
— +  2  7m  V7\l  +  7'  ~  ~ 

v7m  y  V  2?7m  ^ +  27, 


1  f  1 

aa  VoV? 

,  —  +  27m 


1 — LA  1Jhi_rXiV{ 

4fjm  A  +  2rjm  J 


A  .  ^ 

,  —  +  27, 

1 7m 


2 _ 1  ^-2^  „ 

67m  ^  +  27m2 


Plugging  this  approximation  back  into  our  expression  for  the  integral  yields: 


drj 


f)m-£ 


1-7 2  -  Ain 
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7« 
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In  turn,  substituting  this  integral  into  the  expression  for  the  deflection  angle  yields: 


%  =  7r-2A  (p 


r/m 

=  n  -  2sin_1  7  -  2  | 

<7« 


dr/ 


I-72  -Ain— 

I  7oo 


2^ 2 
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K7m 
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1  ?lA^lL£ 

67 m  A  +  27m2 


(4.36) 


(4.37) 
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4.3  Reference  Tables 


Since  we  cannot  invert  the  previous  integral  to  obtain  an  impact  parameter  as  a  direct 
function  of  a  given  deflection  angle,  we  produce  reference  tables  containing  values  of  % 
for  a  range  of  values  of  both  b  and  X.  During  our  calculations,  given  a  value  of  both  b 
and  X,  we  can  interpolate  the  corresponding  value  of  %. 

4.3.1  Positive  Tether 

For  very  low  values  of  X,  the  positive  tether’s  sheath  produces  deflection  angles  similar 
to  those  predicted  by  the  soft-angle  approximation.  For  extremely  low  values  of  X, 
however,  the  total  value  of  the  deflection  angle  is  dwarfed  by  the  error  in  the  integral 
approximation,  resulting  in  negative  values,  though  of  negligible  order. 

4.3. 1.1  Positive  Tether  -  Minimum  Radius 


Minimum  Radius  for  Ions  Entering  Positive  Tether 


lambda 


Figure:  4-4:  Minimum  Radius,  Positive  Tether 

For  very  low  values  of  X,  the  minimum  distance  approaches  zero,  as  the  incoming 
particles  are  deflected  much  more  rapidly. 


96 


4.3. 1.2  Positive  Tether  -  Deflection  Angle 


Deflection  Angle  for  Positive  T ether 


2  4  6  8  10 


lambda 

Figure:  4-5:  Deflection  Angle,  Positive  Tether 

For  any  given  value  of  lambda,  the  smaller  the  impact  parameter,  the  greater  the 
deflection  angle.  This  makes  sense,  since  we  expect  the  tether  to  produce  a  greater  effect 
upon  the  ion  when  it  makes  a  closer  pass.  For  the  further  portions  of  the  graph 

4.3. 1.3  Positive  Tether  -  Deflection  Angle  Focused 

The  problem  with  this  graph  is  that  the  high  gradient  of  the  deflection  angle  for  low 
values  of  X  makes  interpolation  imprecise  for  low  values  of  b  and  A.  Our  solution  is  to 
adjust  our  sampling  positions  such  that  it  is  focused  about  the  area  of  highest  change.  For 
b,  we  space  our  sampling  points  inverse-exponentially  from  b=0: 
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b  (fraction  of  radius  sheath) 


1 


Deflection  Angle  for  Positive  Tether,  Focused 


exp(-l) 

exp(-2) 

exp(-3) 

exp(-4) 

exp(-5_) 

exp(-6) 


lambda 


Figure:  4-6:  Deflection  Angle,  Positive  Tether,  Focused 
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4.3.2  Negative  Tether 


4.3.2. 1  Negative  Tether  -  Minimum  Radius 

First  off,  we  look  at  the  values  of  rm: 


Mimimum  Radius  for  Ions  Entering  Negative  Tether  Sheath  at  2000km 


123456789  10 


lambda 

Figure  4-7:  Minimum  Radius,  Negative  Tether 

Naturally,  since  this  is  an  attracting  tether,  the  minimum  radii  will  be  much  smaller  than 
those  for  the  positive  tether.  For  very  low  values  of  X  (i.e.  for  ions  with  very  high 
energy),  the  tether  has  little  chance  to  exert  a  force  on  the  ion,  and  thus  the  minimum 
radius  approaches  the  initial  impact  parameter  as  X  approaches  zero. 

4.3.2.2  Negative  Tether  -  Deflection  Angle 

When  we  attempt  to  replicate  the  data  with  the  negative  tether,  we  initially  expect  our 
results  to  be  vaguely  similar  to  our  results  to  our  positive  tether,  with  smaller  impact 
parameters  resulting  in  deflection  angles  of  larger  magnitude,  keeping  in  mind  that  the 
deflection  angles  will  now  be  measured  as  negatives.  However,  when  we  calculate  the 
deflection  angle,  we  discover  an  unusual  phenomenon:  for  any  lambda,  the  deflection 
angle  as  a  function  of  impact  parameter  is  no  longer  monotonic. 
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Deflection  Angle  for  Negative  Tether 


2  4  6  8  10 

lambda 


Figure  4-8 :  Deflection  Angle,  Negative  Tether 

For  the  lower  values  of  X,  the  deflection  angle  once  again  appears  to  obey  the  soft-body 
approximation  and  approach  zero,  yet  it  is  apparent  from  the  contour  graph  that  the 
deflection  angle  magnitude  first  increases  with  increasing  impact  parameter  and  after 
some  turning  point  starts  decreasing  again.  This  turning  point  appears  to  increase  with  X. 
The  scale  of  the  graph  makes  the  contour  lines  for  high  values  of  b  difficult  to  read,  but 
we  know  that  for  any  finite  value  of  lambda,  the  deflection  angle  must  equal  zero  if  the 
impact  parameter  equals  the  radius  of  the  sheath.  Assuming  there  are  no  discontinuities 
in  the  graph,  that  must  mean  that  for  every  value  of  X,  there  exists  a  value  bmax  such  that 
the  magnitude  of  the  deflection  angle  is  maximum  for  that  value  of  X.  For  low  values  of 
X,  this  value  is  a  fairly  small  fraction  of  the  sheath  radius.  For  larger  values,  the  distance 
between  this  value  and  the  sheath  radius  is  tiny. 

4.3.2.3  Negative  Tether  -  Deflection  Angle  Focused 

For  computational  purposes,  we  take  an  exponential-scale  sample  of  the  data  to  zoom  in 
on  the  impact  parameters  closest  to  the  radius  sheath. 
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Deflection  Angle  for  Negative  Tether,  Focused 

1-exp(-5) 


1-exp(-4) 


1-exp(-3) 


1-exp(-2) 


l-exp(-l) 
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lambda 

Figure  4-9:  Deflection  Angle,  Negative  Tether,  Focused 

The  exponential  close-up  verifies  our  previous  analysis,  as  we  can  now  more  clearly 
define  the  values  of  bmax  for  the  lower  values  of  X.  For  high  values  of  X,  even  though  we 
know  bmax  exists,  it  is  even  greater  than  our  sampling  increment  closest  to  the  sheath 
radius.  When  we  compute  the  ranges  for  the  tether,  the  range  of  b  greater  than  bmax 
becomes  negligible. 

In  order  to  determine  the  source  of  this  odd  behavior,  we  compare  the  above  graphs  of 
minimum  radius  and  deflection  angle.  For  a  given  X,  there  exists  some  value  bmax  such 
that  the  maximum  possible  deflection  angle  is  achieved.  For  biow<  bmax,  the  ion  makes  a 
closer  pass  to  the  tether,  but  the  pass  doesn’t  last  as  long.  What  appears  to  be  happening 
is  that  because  the  ion  has  a  lower  value  of  b,  more  of  the  force  exerted  upon  it  by  the 
tether  is  directed  in  the  direction  parallel  to  its  velocity  until  it  passes  very  close  to  the 
tether.  For  bmax,  the  tether  potential  exerts  more  force  in  the  radial  direction,  translating 
into  a  weaker  but  longer-lasting  centripetal  force. 

If  you  could  position  an  ion  directly  into  the  sheath  such  that  its  velocity  vector  was 
exactly  perpendicular  to  the  radial  direction,  there  would  exist  a  value  X  such  that  the 
force  on  the  ion  would  be  precisely  equal  to  the  centripetal  force  required  for  a  stable 
orbit,  and  the  ion’s  deflection  angle  would  be  infinite.  We  know  this  to  be  an  impossible 
scenario  for  the  ions  in  question  because  the  only  way  an  incoming  ion  can  enter  the 
sheath  at  precisely  90°  would  be  if  the  impact  parameter  equaled  the  sheath  radius,  at 


_Q 
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which  point  no  deflection  occurs  at  all.  However,  such  a  trajectory  can  be  partly 
mimicked  if  the  ion  approaches  the  sheath  at  an  angle  and  speed  that  maximizes  angular 
acceleration  while  minimizing  linear  acceleration. 


Figure  4-10:  Maximum  Impact  Parameter 


With  this  in  mind,  let  us  examine  the  integral  for  the  deflection  angle  itself.  Of  particular 
note  is  the  integral  defining  the  angle  transversed  by  the  ion  while  within  the  sheath: 


A<p  =  $ 
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1  — — r  — /tin  — 
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(4.38) 


The  larger  this  integral  becomes,  the  greater  the  magnitude  of  the  deflection  angle 
becomes.  Thus,  let  us  examine  how  increases  in  b  and  A  affect  the  size  of  this  integral. 

The  limits  of  the  integral  are  rm  and  oo.  Set  A,  as  a  constant.  We  already  know  that  as  b 
increases,  rm  increases,  and  thus  the  range  of  integration  decreases.  However,  if  we  look 
at  the  content  of  the  integral,  we  see  that  an  increase  in  b  causes  an  increase  in  the 
integral’s  interior.  Our  adjustments  in  b  thus  produce  two  counteracting  effects  on  the 
integral.  Below  bmax,  the  increasing  effects  dominate,  and  the  integral  increases  with  b. 
Above  bmax,  the  limit  restriction  dominates,  and  the  integral  decreases  with  b. 
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Chapter  5 

Results  and  Discussion 


5.1  Tether  Scattering  Flux 

For  each  of  the  tether  altitudes,  we  have  been  able  to  calculate  not  only  the  total 
scattering  flux  per  unit  length  of  tether,  but  also  the  total  influx  of  particles  into  the 
sheath  per  unit  length.  Thus  we  are  able  to  compute  not  only  the  scattering  flux  of  the 
entire  tether,  but  also  the  efficiency  of  the  tether.  Note  once  again  that  we  will  be 
focusing  specifically  on  the  population  of  high-energy  ions,  as  low-energy  ions  are  of 
little  concern  to  us,  and  our  calculations  cannot  be  used  to  determine  electron  fluxes 
accurately. 

The  first  term  we  will  evaluate  will  be  the  total  flux  of  particles  entering  the  sheath, 
which  we  shall  call  the  sheath  influx.  We  take  our  particle  influx  equation  from  4.12  and 
obtain: 
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(5.1) 


Since  we  restricting  our  focus  to  only  the  population  of  high-energy  particles,  we  set  the 
temperature  T  to  correspond  to  particles  of  the  energy  IMeV,  such  that  T  =  106  •  11600. 
The  ambient  density  term  in  the  distribution  function  could  easily  be  extracted  from  the 
flux  integral  before  calculation: 
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(5.2) 


Since  separating  the  term  proves  useful  for  our  analysis  later,  we  display  the  particle 
influx  as  a  product  of  the  ambient  density  and  the  remaining  integral,  which  is  now  just 
the  volume  influx,  or  y: 
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Substituting  the  sheath  radii  for  each  altitude  at  each  tether  potential  yields  the  following 
volume  influxes: 


Altitude  (km) 

2000 

4000 

6000 

8000 


Volume  Influx  (mV1) 
Positive  Tether 
1.5153x10" 
1.4678x10" 
1.3792x10" 
1.3181x10" 


Table  5-1:  Volume  Influx ,  Positive  and  Negative  Tethers 


Volume  Influx(m3s_1) 
Negative  Tether 
2.9786xl013 
2.8819xl013 
2.7016xl013 
2.5775xl013 


As  we  increase  the  tether’s  altitude,  the  loss  cone  angle  decreases,  the  effects  of  which 
are  twofold.  First,  since  we  assume  the  velocities  of  all  particles  in  the  magnetosphere  lie 
outside  the  loss  cone,  decreasing  the  loss  cone  angle  increases  the  range  of  radial  velocity 
g  for  the  incoming  ions.  Since  the  flux  incorporates  an  extra  g  term  into  the  normalized 
distribution  integral,  the  increased  range  of  g  makes  for  a  greater  range  of  angles  for  the 
ions  to  penetrate  the  sheath,  thus  increasing  the  flux.  However,  this  effect  is  counteracted 
by  the  fact  that  the  tether  has  to  redirect  incoming  ions  into  a  much  smaller  loss-cone, 
which  it  misses  more  frequently.  That  is,  even  though  the  total  number  of  incoming 
particles  increases  with  altitude,  the  total  number  that  is  actually  scattered  into  the  loss 
cone  decreases. 


To  calculate  the  particle  scattering  flux,  we  similarly  split  up  the  integral  into  a  product  of 
the  ambient  density  and  the  volume  scattering  flux: 

r scatter  =  \\\  g2  '  f Lc{V  z  >  g  >9)' Wb{V  z  ’  g’9)' LT  ’ dgdOdv ZV 

Outside 
Loss  Cone 


V1V-  1 


Outside 
Loss  Cone 


COS  cc 


m 


L  V 


2nkT 


e  2kT  ■  w. 


( vz ,  g,  0)  ■  Lt  ■  dgd 9dvzm 


3-1 


f  ’  =  ri  y 

scatter  oo/  scatter 


(5.4) 
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From  this  we  derive  our  volume  scattering  fluxes: 


Altitude  (km) 

Volume  Scattering  Flux 

Volume  Scattering  Flux 

Positive  Tether  (m  s'  ) 

Negative  Tether  (m  s' ) 

2000 

8.6978xl09 

1.8130xl012 

4000 

3.9829xl09 

8.6175x10" 

6000 

1.3792x10" 

4.4095x10" 

8000 

1.3181x10" 

1.5059xl010 

Table  5-2:  Volume  Scattering  Flux 

We  can  also  calculate  the  tether  “efficiencies”  by  dividing  the  volume  influx  by  the 
volume  scattering  flux  for  each  altitude  at  each  tether: 

Altitude  (km) 

Efficiency  (%) 

Efficiency  (%) 

Positive  Tether 

Negative  Tether 

2000 

5.74 

6.09 

4000 

2.71 

2.99 

6000 

1.45 

1.63 

8000 

0.02 

0.06 

Table  5-3:  Tether  Efficiency 

The  most  obvious  result  from  our  analysis  is  the  rate  at  which  tether  efficiency  is  reduced 
as  we  increase  altitude. 


5.2  Remediation  Time 

To  determine  how  long  the  tether  would  take  to  deplete  a  certain  region  of  space,  we  first 
observe  that  its  scattering  flux  is  directly  proportional  to  the  ambient  plasma  density. 

Assume  we  wish  to  thin  a  certain  region  of  the  magnetosphere  by  a  factor  of  (1  (that  is, 
reducing  the  magnetosphere’s  population  to  one-tenth  its  original  value  would  translate  to 
P=10).  Say  that  the  tether  has  a  scattering  flux  of  yn.  Within  a  single  unit  of  time,  the 
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tether  thus  scatters  yn  ions.  Now  suppose  we  isolate  the  tether  within  a  given  volume  V, 
which  initially  contains  a  total  number  of  particles  N,  such  that  the  ambient  density  is  n. 
To  calculate  the  depletion  rate,  we  first  calculate  the  factor  by  which  the  ambient  density 
is  reduced  over  a  unit  of  time  At.  We  start  by  defining  the  density  before  and  after  this 
unit  of  time: 


n,  = 


N_ 

y 


n/  = 


(N-W scatter  At) 


(5.5) 


Now  obtain  An  by  subtracting  the  initial  density  from  the  final  density: 


An  =  nf  - 


(N-ny^At)  N 
V  V 


(5.6) 


Divide  by  the  given  unit  of  time  to  obtain  a  differential  equation  defining  the  density: 


dn  _  An  _  y  ■  n  (5.7) 

dt  At  V 


Integrate  both  sides  to  obtain 


n  =  n 


o 


(5.8) 


Thus,  we  can  define  a  target  density  we  wish  to  achieve  for  this  region,  and  the  following 
tells  us  how  long  this  goal  will  take: 


^ final 


n 


0 


t  = 


\ 

W0 

yn final  y 


(5.9) 


Now  to  define  the  region  we  wish  to  isolate.  Since  we  are  attempting  to  deplete  the 
radiation  belts  via  a  tether  traveling  in  the  equatorial  plane  at  a  certain  altitude,  we  limit 
our  space  to  those  magnetic  field  lines  which  intersect  the  equatorial  plane  near  the 
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altitude  of  our  satellite.  We  approximate  the  magnetic  field  lines  of  a  given  strength  and 
altitude  to  form  a  torus  around  the  Earth. 

The  volume  of  a  torus  is  given  by: 

F„=2  x2Rr2  (5.10) 


R  is  the  radial  distance  from  the  torus  center  to  its  circular  axis,  and  r  is  the  radius  about 
the  circular  axis.  To  approximate  the  shape  of  the  magnetic  field  lines,  we  take  both  R 
and  r  equal  to  half  the  radius  to  the  magnetic  field  line  in  question,  or: 


V, . =  -r2Rj 


eq 


(5.11) 


If  we  want  to  determine  the  volume  encompassed  by  the  magnetic  field  lines  that 
intersect  the  equatorial  plane  within  a  certain  radial  distance  form  our  tether’s  orbit,  we 
calculate  the  volume  of  the  torus  with  the  larger  radius  and  subtract  the  volume  of  that  of 
the  smaller  radius. 


5  —  R  5 

eq, outer  veq,  inner 


(5.12) 


When  we  substitute  (5.8)  into  our  equation  for  mission  time,  we  discover  that  density  by 
itself  plays  no  role  in  the  total  mission  time,  while  the  desired  fraction  depletion  does. 


1  2  3 
^  ^  (^eq, outer 


R . 


eq, inner 


>08) 


(5.13) 


For  a  given  altitude  of  tether  orbit,  let  us  define  our  target  volume  by  the  area  covered  by 
the  magnetic  field  lines  which  intersect  the  equatorial  plane  at  altitudes  within  a  distance 
D  of  the  tether  orbit.  Or: 


R  _  d  _i _  d 

^eq, outer  ^ orbit  '  ^ 

R  =  R  —  D 

Veq, inner  orbit 


(5.14) 


Substitute  this  into  our  time  equation  to  get: 
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1  =±x\(RoMt+Df  -(R.M,-DfMp) 


(5.15) 


Furthermore,  the  satellite  possesses  both  positive  and  negative  tethers,  each  of  cumulative 
length  10  km.  Thus,  the  term  for  the  scattering  rate  is  actually  the  sum  of  the  rates  from 
each  tether,  or: 


=  4  {r^ry%^ADi-(KM,-D)'MP) 


(5.16) 


Let  us  suppose  we  restrict  our  toroidal  space  such  that  it  covers  only  those  magnetic  field 
lines  whose  equatorial  altitudes  are  directly  intersected  by  the  tether’s  orbital  path;  that  is, 
we  take  D  to  equal  10km.  If  we  wish  to  reduce  the  density  within  this  isolated  range  to 
one-tenth  of  its  original  value  (i.e.  set  (3=10),  we  get  the  following  times  at  various 
altitudes: 


Altitude 

T (sec) 

T  (years) 

2000 

1.3135xl07 

0.4165 

4000 

6.0330xl07 

1.9130 

6000 

2.0643xl08 

6.5458 

8000 

1.5722xl010 

498.54 

Table  5-4:  Mission  Time:  D  =  10km 

At  low  altitudes,  the  mission  time  falls  to  within  approximately  one  year,  so  depleting 
this  section  of  the  magnetosphere  seems  a  plausible  task.  At  high  altitudes,  however,  the 
mission  time  increases  to  outrageous  proportions,  so  we  might  be  better  off  setting  out 
sights  a  bit  lower. 

Of  course,  these  results  are  unrealistic  in  that  we  cannot  cordon  off  the  area  of  space 
surrounding  those  we  have  designated  to  be  part  of  V  without  particles  from  the 
surrounding  areas  randomly  scattering  into  it  and  repopulating  it.  Thus,  let  us  expand  our 
borders  by  a  factor  of  10  and  define  D  =  100km.  Our  results  are  now: 


Altitude 

T  (sec) 

T (years) 

2000 

1.3135xl08 

4.1652 

4000 

6.0331xl08 

19.131 
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6000 


65.459 


2.0643  xlO9 

8000  1.5722x10"  4985.5 

Table  5-5:  Mission  Time:  D  =  100km 
Expanding  it  to  D  =  1000km: 


Altitude 

T (sec) 

T  (years) 

2000 

1.3140xl09 

41.847 

4000 

6.0461  xlO9 

191.72 

6000 

2.0669xl010 

655.39 

8000 

1.5730xl012 

49878. 

Table  5-6:  Mission  Time:  D  =  1000km 

The  volume  increases  roughly  linearly  as  a  result  of  our  range  expansion,  so  naturally  the 
depletion  time  is  also  linearly  increased  as  a  result. 
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Chapter  6 

Conclusion  and  Recommendations 


6.1  Conclusion 

The  prospects  of  a  single  tether  satellite  depleting  a  substantial  portion  of  the  radiation 
belts  now  is  temporally  feasible  for  low  altitudes  around  2000km,  but  not  for  higher 
altitudes  around  8000km.  When  designing  missions  for  magnetospheric  remediation  of 
high-energy  ions,  a  single  tether  should  suffice  for  the  low  altitudes,  but  an  array  of 
several  dozen  satellites  might  be  required  for  the  higher  altitudes. 

However,  even  for  the  lower  altitudes,  the  tether  satellite  can  sustain  depletion  only  for 
short  spans  of  altitude.  To  deplete  a  range  of  altitude  on  the  order  of  1000km,  one  would 
need  an  array  of  satellites  constantly  functioning  at  varying  altitudes.  Note  that  we  have 
yet  to  consider  the  effects  of  natural  replenishment  from  cosmic  neutrons,  so  our 
estimates  would  appear  to  be  best-case  scenarios  which  must  be  tempered  with  future 
research. 

Regarding  the  implausibility  of  remediation  missions  at  higher  altitudes,  on  of  the 
greatest  hindrance  is  the  loss  cone  angle,  which  gets  terribly  narrower  as  we  increase 
altitude,  resulting  in  reduced  efficiency,  especially  above  6000km.  At  the  higher 
altitudes,  even  if  every  particle  to  interact  with  the  tether  satellite  were  to  be  deflected 
into  the  loss  cone,  the  number  of  particles  to  be  deflected  would  still  require  at  least 
decades  for  the  remediation  requirements  missions  we  wish  to  fulfill.  Also,  while  our 
calculations  do  not  consider  possible  influx  rate  increases  as  a  result  of  our  tether’s 
orbital  velocity,  yet  any  increase  of  influx  on  the  ram  end  would  likely  be  offset  at  least 
in  part  by  a  wake  on  the  tail  end,  so  any  net  gain  in  scattering  rate  would  not  in  itself 
deem  feasible  any  high-altitude  missions. 

One  method  of  rendering  high-altitude  tethers  more  efficient  would  be  to  increase  the 
sheath  size,  which  would  necessitate  increasing  either  the  sheath  radius  or  the  sheath 
length.  Increasing  the  latter  is  a  simple  matter  of  increasing  the  overall  length  of  the 
tether,  but  the  resulting  increase  in  ion  influx  would  be  linear.  If  we  wanted  to  increase 
the  ion  influx  by  two  orders  of  magnitude,  we  would  increase  the  tether  length  similarly, 
from  a  10-km  tether  to  a  1000-km  tether.  Since  the  tether  radius  term  is  logarithmic 
while  the  tether  potential  term  is  not,  changes  in  the  latter  overshadow  changes  in  the 
former.  Furthermore,  the  added  tether  length  would  also  result  in  increases  in  charge 
collection,  which  would  increase  the  deforming  Lorentz  forces  and  reduce  the  tether 
radius. 
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If  we  choose  to  increase  the  sheath  radius,  we  still  encounter  problems  when  considering 
the  power  and  structural  limitations.  Increasing  the  tether  potential  by  several  orders  of 
magnitude  would  easily  increase  the  sheath  radius  to  a  more  acceptable  size,  but  the 
tether  radius  would  have  to  be  reduced  by  an  even  greater  factor  to  compensate  for  the 
power  limitation: 


rT 


^max  /  m, 

na>L]j-2  (e<V)3 


(6.1) 


Increasing  the  tether  potential  by  a  factor  of  100  thus  decreases  the  tether  radius  by  a 
factor  of  1000,  making  our  tether  radius  roughly  0.5pm.  Our  design  for  a  0.5mm  tether 
was  tenuous  enough  as  was,  even  if  our  design  employed  tungsten  steel,  but  if  our  tether 
is  meant  to  connect  two  100-kg  power  supplies,  0.5-pm  is  much,  much  more  likely  to 
snap.  We  could  increase  the  tether  radius  by  also  decreasing  the  length  of  the  tether,  but 
that  would  linearly  reduce  the  size  of  the  sheath,  and  the  ion  influx,  thus  defeating  our 
intended  purpose. 


6.2  Recommendations 

For  future  research,  we  would  recommend  several  avenues  for  increasing  the  efficiency 
of  scattering  tether  mission  designs.  If  one  were  to  employ  the  series  tether  design 
outlined  in  this  thesis  and  expect  a  feasible  mission  plan,  one  would  be  required  to 
analyze  different  models  for  tether  arrays  in  various  orbits  and  maximize  their  cumulative 
depletion  rates.  One  should  also  analyze  and  compare  various  other  designs,  including 
tether-clusters  and  tape  tethers,  to  determine  the  most  efficient  and  effective  scattering 
methods.  Further  research  is  also  required  to  determine  the  effects  of  each  of  these  tether 
designs  upon  the  natural  replenishment  rates  of  the  radiation  belts. 

Whichever  tether  design  one  should  settle  upon,  one  must  also  determine  the  feasibility 
of  its  design  with  regards  to  structural  mechanics.  Our  analysis  of  Lorentz  forces  upon 
the  series  tether  concludes  that  non-negligible  forces  are  exerted  lengthwise  on  either  end 
of  the  tether,  and  whatever  material  is  used  to  construct  the  satellite  must  be  able  to  hold 
its  shape  under  these  forces  for  consistent  scattering  operation.  We  calculated  our 
scattering  rates  on  the  assumption  that  a  0.5-mm  tether  could  survive  collisions  with 
heavy  ions  and  cosmic  debris,  but  this  assumption  must  be  verified  before  we  can  cast 
away  the  concern  that  the  tether  may  snap  too  easily  once  deployed.  We  also  assumed 
that  attaching  two  100-kW  power  supplies  to  the  series  tether  was  a  plausible  design,  yet 
this  electronic  configuration  must  too  be  verified  for  feasibility  over  long  periods  of  time. 


Ill 
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APPENDIX  3:  PITCH  DIFFUSION  BY  RESONANT  INTERACTION 

WITH  WHISTLERS 

Radiation  Belt  Remediation  methods 

In  order  to  reduce  the  density  of  energetic  particles  inside  the  radiation  belts  we  have  to 
make  them  enter  their  loss  cone.  There  are  three  ways  of  doing  so,  or  any  combination  of 
them: 

•  To  increase  their  parallel  energy  /  momentum. 

•  To  reduce  their  perpendicular  energy  /  momentum. 

•  To  reduce  their  incidence,  rotating  their  velocity  without  any  change  in  their 
energy. 

The  use  of  space  tethers  for  modifying  the  particles’  trajectories  has  already  been 
considered.  We  can  distinguish  between  two  main  approaches. 

The  first  is  using  high  voltage  electrostatic  tethers.  The  electric  potential  around  the  tether 
deviates  the  particles,  making  some  of  them  enter  their  loss  cone.  The  tether’s  orbital 
motion  allows  it  to  cover  a  significant  volume  inside  the  magnetosphere.  A  preliminary 
study  for  such  a  system  has  already  been  done  [ Christopher  F.  Zeineh ],  and  it  shows  how 
the  characteristic  time  for  remediation  is  of  the  order  of  the  year  for  a  10  kW  system  at  an 
altitude  of  L=2.  For  higher  values  of  L  the  effectiveness  of  the  system  decreases  quite 
dramatically. 

The  second,  which  is  the  main  objective  of  this  preliminary  study,  is  to  use  space  tethers 
as  VLF/ELF  antennas.  These  antennas  emit  whistler  waves  which  can  perturb  the  motion 
of  the  energetic  trapped  particles,  creating  a  net  diffusion  towards  their  loss  cone  and  thus 
removing  them  from  the  radiation  belts.  We  expect  the  characteristic  remediation  times  to 
be  significantly  smaller,  since  such  a  system  does  not  act  at  a  local  but  at  a  global  scale. 
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Theoretical  Background 

In  this  chapter  we  will  analyze  the  wave-particle  interaction  process  and  explain  what 
will  be  the  theoretical  frame  used  to  solve  the  problem  and  obtain  the  characteristic 
precipitation  time. 

Electron  cyclotron  waves 


The  first  idea  that  comes  to  our  mind  to  perturb  the  energetic  particles’  motion  is  to  emit 
waves  at  the  resonant  frequency  of  their  Larmor  motion.  In  doing  so,  the  electric  and 
magnetic  fields  seen  by  the  particles  will  be  quasi-static  and  there  will  be  a  significant 
energy  exchange  between  both  of  them.  But,  since  this  frequency  is  independent  of  the 
particle’s  energy,  this  is  not  a  good  idea.  Using  these  waves,  the  background  plasma,  with 
a  density  several  orders  of  magnitude  higher  and  energies  of  the  order  of  0.4  eV,  will 
absorb  them,  and  such  a  system  will  do  nothing  but  heating  this  plasma. 

To  verify  the  importance  of  this  dumping,  we’ve  taken  this  dispersion  relation,  given  by 
the  kinetic  theory  of  plasma  perturbations 


n 2 =  1  +  - 
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and,  being  the  background  plasma  non  relativistic,  we  can  consider  their  mean  velocity  as 


v„  = 


2  kBTe 
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(2.2) 


The  plasma  dispersion  function  can  be  approached  by  the  following  development, 
because  we  are  near  the  cyclotron  frequency 
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and  the  dispersion  relation  can  thus  be  written 
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To  obtain  an  estimate  of  the  penetration  distance  of  the  waves  inside  the  background 
plasma,  we  consider  complex  solutions  for  the  equation  (2.4).  The  inverse  of  the 
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imaginary  part  of  the  wave  number  corresponds  to  a  characteristic  penetration  length. 
Solving  equation  (2.4)  is  equivalent  to  obtaining  the  roots  of  the  following  complex 
polynomial. 


E,4  —  v2E,2 

where  the  variables  are 


^co  V'  ^ 


pe 


z'  \ 


v<;  + 


CO 


pe 


f  r.\ 


KVeJ 


V 


(v-l)=0 


(2.5) 
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Obtaining  the  roots  for  L=2  and  frequencies  near  the  cyclotron  one,  we’ve  plot  the  results 
and  we  see  that  after  some  tens  of  meters,  the  waves  are  absorbed  by  the  background 
plasma. 


Space  Constant 
L=2 


Figure  2-1.  Characteristic  penetration  distance  for  an  electron  cyclotron  wave  as  a 

function  of  its  frequency  for  L=2 

The  use  of  such  waves  is  of  no  interest  at  all,  because  all  the  power  is  used  to  heat  the 
background  plasma,  instead  of  perturbing  the  motion  of  the  more  energetic  particles. 
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VLF/ELF  emissions 


Low  frequency  modes,  like  the  whistlers,  can  induce  particles’  precipitation.  In  nature, 
these  modes  can  be  excited  by  lightning,  but  we  can  emit  such  waves  with  antennas  of 
the  proper  length.  The  influence  of  this  and  other  low  frequency  modes  on  the 
magnetospheric  particles’  distribution,  and  its  possible  use  for  radiation  belt  remediation, 
has  already  been  considered  in  several  studies  like  [Abel  &  Thorne,  1998 ]  and  [Inan, 
Bell,  Bortnik  &  Albert,  2003],  We  will  see  now  how  to  calculate  the  effects  of  such 
waves  on  the  energetic  populations  of  the  ionosphere. 

Cyclotron  resonance 


The  wave  particle  interaction  at  low  frequencies  can  be  explained  by  the  existence  of 
Doppler  resonances.  To  obtain  the  resonance  condition,  we  can  start  by  writing  the 
motion  equations  for  an  electron  under  the  influence  of  the  electric  and  magnetic  fields  of 
an  electromagnetic  wave.  The  electrons  considered  are  relativistic  because  their  energies 
are  of  the  order  of  1  MeV.  Once  we  take  into  account  the  magnetic  moment  conservation 

2 

=  cste  (2.7) 

Bo 

and  if  we  take  the  z  axis  as  the  direction  of  the  Earth’s  magnetic  field,  we  obtain 
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here,  the  angle  between  the  perpendicular  momentum  and  the  wave’s  magnetic  field  ,  Bv 
,  is  n  -  op . 


If  the  electrons’  motion  is  to  be  modified  significantly  by  the  wave’s  fields,  their  effect 
must  cumulate  for  a  certain  period  of  time,  and  the  resonance  condition  must  be 

op  =  0  (2.9) 
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so  that  the  electron  does  not  see  rapid  oscillating  fields,  which  mean  effect  is  negligible. 
The  last  term  of  the  third  equation  in  (2.8),  which  gives  the  change  in  the  phase  angle,  is 
negligible  for  the  conditions  given  in  the  Earth’s  magnetosphere.  Therefore 

co =  (2.10) 

mj  y 

and  (2.10)  will  we  used  as  resonance  condition  for  all  our  calculations. 

Cyclotron  waves  whose  frecuency  is  small  in  comparison  with  the  electron 
gyrofrequency,  Qe,  can  thus  resonate  with  the  energetic  particles,  perturbing  their 

motions  in  a  significant  way.  The  particles’  pitch  angle  will  change  stochastically,  some 
times  it  will  increase  and  other  it  will  decrease,  depending  on  the  phase  angle  value  at  the 
moment  of  the  interaction.  Nonetheless,  given  the  fact  that  the  distribution  is  a  hollow 
one,  because  of  the  loss  cone,  there  are  more  particles  with  higher  values  of  pitch  angle, 
and  there  will  be  a  net  diffusion  towards  the  loss  cone. 

Quasilinear  Theory 

The  instabilities’  saturation  process  inside  a  plasma  can  be  described  as  a  continuous 
diffusion  in  the  velocity  space,  in  which  a  zero  order  distribution  function  evolves  slowly. 
The  diffusion  rate  is  proportional  to  the  sum  of  the  squares  of  the  modes  obtained  by  the 
linear  theory,  thus  the  name  of  quasilinear  theory.  When  applying  this  theory  we  are 
making  implicitly  two  assumptions: 

First,  the  amplitudes  of  the  perturbations  must  be  small  enough.  For  small  enough  we 
understand  that  the  results  obtained  for  the  zero  order  distribution  function  given  by  the 
linear  theory  must  be  correct.  We  will  average  this  distribution  function  in  space 

/ofrO=(/ofcv,f)). 

Second,  the  wave  spectrum  must  be  dense  enough  so  that  any  coherence  between  the 
different  modes  is  destroyed  by  phase  mixing  in  the  time  scale  of  the  plasma  parameters 
variation. 

If  we  start  from  the  Vlassov  equation  for  the  plasma 
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(2.11) 


— +  v -V/  +  — V- 
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we  can  average  it  over  a  certain  number  of  periods  in  space  and  time,  and  also  we  can 
average  it  over  the  gyro  motion  around  the  Earth’s  magnetic  field,  B0 .  Doing  so,  we  can 
rewrite  (2.1 1)  under  this  new  form  [. Kennel  &  Engelmann,  1966 ] 
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V2  v± 

And  after  some  math,  we  arrive  to  a  diffusion  equation  for  f0  (E,  L,t,  a0).  Because  of  the 

conditions  given  in  our  case  of  interest,  the  energy  exchange  between  waves  and  particles 
is  negligible,  and  the  diffusion  equation  simplifies  into  a  pure  diffusion  equation  for  pitch 
angle. 


dt 
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sin2aor(a0)  da0 


sin  2aor(a0  )l)aa  (a0 ) 


Ml 

da 


(2.13) 


o  J 


In  (2.13)  we  cannot  see  any  source  of  particles.  That  is  the  case  for  an  artificial  radiation 
belt,  created  for  instance  by  means  of  a  nuclear  detonation  at  high  altitude.  If  we  wanted 
to  study  the  dynamics  of  the  natural  radiation  belts,  the  Van  Allen  belts,  we  would  need 
to  model  the  particles’  sources,  but  this  is  out  of  the  scope  of  this  study.  But  for  the  case 
of  the  inner  Van  Allen  belt,  equation  (2.13)  could  be  taken  as  an  approximation,  for  the 
dynamics  of  its  sources  are  very  slow. 


To  solve  (2.13),  the  easiest  way  is  to  assume  that  we  can  separate  both  variables,  time 
and  pitch  angle 

f0(t,a0)=N{t)g(a0)  (2.14) 

N  represents  the  equatorial  electron  density,  for  a  given  value  of  L  and  E,  and  g  gives  us 
the  shape  of  the  distribution  in  pitch  angle.  We  consider  the  density  to  be  constant  along 
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the  field  lines,  and  g  constant  along  the  time  evolution  of  the  distribution,  or  at  least  that 
the  shape  of  the  distribution  changes  very  slowly  in  comparison  with  the  density 


1  clNjt  )  _ _ 1 _ d_ 

N(t )  dt  sin2a07’(a0)g(a0)  da0 


sin  2a07’(a0  )Daa  (a0 ) 


Jg(a0) 


da 


(2-15) 


o  J 


here  x  is  the  characteristic  time  scale  for  the  equation,  at  it  is  representative  of  the  time 
needed  to  clean  the  L  shell  in  which  we  find  our  antenna 


f(a0 )  is  an  approximate  function  which  gives  us  the  bounce  period  as  a  function  of  the 
equatorial  pitch  angle,  taking  the  Earth’s  magnetic  field  as  dipolar. 

r(a0)  =  1.3802 -0.3 198(sina0  +^/sina0 )  (2.16) 

Several  authors  [Lyons  &  Thorne,  1972],  [Abel  &  Thorne,  1998],  solve  (2.15)  using 
iterative  methods.  Nevertheless,  these  methods  can  give  some  convergence  problems  if 
the  diffusion  coefficient  approaches  zero  for  certain  pitch  angle  values. 

That  is  why  we  have  decided  to  rewrite  the  problem  under  the  form  of  an  eigenvalue 
problem. 
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(2.17) 


The  operator’s  matrix  dimension  is  equal  to  the  number  of  pitch  angles  chosen  to 
discretize  the  electrons’  distribution  function.  In  principle,  the  negative  eigen  value  which 
yields  the  shortest  precipitation  time  is  the  one  we  are  looking  for,  and  its  corresponding 
eigen  vector,  the  distribution  in  pitch  angle,  g. 


But  numerical  instabilities  can  make  more  difficult  the  identification  of  the  eigen  value 
we  are  looking  for.  In  this  case,  easily  identifiable  if  the  corresponding  g  does  not  satisfy 
the  conditions,  we  have  to  verify  each  one  of  the  eigen  values,  in  a  process  that  can  be 
very  time  consuming. 
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We  have  to  impose  two  boundary  conditions  to  g.  The  first  one  is  that  the  distribution 
function  must  be  zero  at  the  loss  cone.  The  second  one  is  that  the  distribution  function 

must  be  symmetrical  with  respect  to  ^ 


g(<* 


0  ,LC 


dg 


da. 


)-0 

=  0 


(2.18) 
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Diffusion  Coefficients 


To  solve  (2.16)  we  need  to  know  the  diffusion  coefficient  as  a  function  of  the  equatorial 
pitch  angle.  This  coefficient  will  be  the  sum  of  the  coefficients  for  each  resonance. 

Dja„)=  I»0)  r  =  0,±1,±2,...  (2.19) 

r 

Here  r  =  0  corresponds  to  the  Landau  resonance.  The  value  of  the  bounce  averaged 
Landau  coefficient  can  be  obtained  trough 

D, K)=  ,  r  i  f"  cosa  •  sin4  a  -  cos7  X  dX  (2.20) 

7\a0)cos  a0  J0  2  B  y  to 

A  is  a  normalization  constant  that  can  be  calculated  if  we  take  the  energy  of  the  wave  to 
be  distributed  in  frequency  and  propagating  angle  following  a  Gaussian. 

( co-to„  V  T  tan  9-tan  0m  f 

B2w  OC  e  ^  8(0  )  -e  ^  tan8e  )  (2.21) 

We  call  0  the  angle  between  the  wave  vector,  k  ,  and  the  local  magnetic  field.  The  rest 
of  the  parameters  needed  to  determine  the  Gaussians  are  defined  by  the  emitter 
properties. 

The  value  for  this  constant  is 
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where  we  have  introduced  a  new  variable,  x  =  tan  0 ,  and  an  low  frequency  cut  off  at 
oo, 


JLC  ■ 


W0  is  the  weighting  function  for  the  Landau  resonance. 
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vF(0)=  exp 

®o(e)= 


x-  tan0fl 
V  tan  80 

co 
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(2.24) 
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and  Jt  are  Bessel  functions  of  the  first  kind,  whose  argument  is 
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The  integration  limit  is 
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and  this  limit  allows  us  to  avoid  all  the  frequencies  under  de  cut  off.  All  these  expressions 
depend  on  two  different  linear  momenta.  The  first  is  the  particles’  linear  momentum 
parallel  to  the  local  magnetic  field 

i^Pcosa  (2.26) 

The  second  is  the  parallel  momentum  required  for  the  first  order  resonance,  m .  We  can 
obtain  the  latter  trough  the  general  resonance  condition 

rCl, 
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For  the  first  order  resonance,  r  =  -1 
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Finally,  all  the  expressions  that  must  be  integrated  in  Wo  should  be  evaluated  for  a 
frequency  obtained  from  the  Landau  resonance  condition 
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For  the  rest  of  the  resonances,  the  bounce  averaged  diffusion  coefficient  can  be  obtained 
as  follows 
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<Me)  =  ^  [(1  +  cos  0)  •  Jr+1  +  (1  -  cos  0)  •  Jr  l  f  (2.32) 

the  argument  of  the  Bessel  functions  is  now  [r  tan  ax],  and  to  avoid  taking  into  account 
the  effect  of  frequencies  below  the  cutoff,  we  have  to  limit  the  integration  from 


(2.33) 


Finally,  all  the  parameters  involved  in  the  integration  of  the  weighting  function  Wr  must 
be  evaluated  for  the  corresponding  resonant  frequency. 
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Modeling 

In  this  chapter  we  will  analyze  the  different  models  used  to  implement  the  theoretical 
expression  given  in  the  preceding  chapter  in  a  computational  code. 

Environment 

Earth’s  Magnetic  Field 

For  the  sake  of  simplicity  we  have  considered  the  Earth’s  magnetic  field  in  the  radiation 
belt  region  as  that  of  a  centered  magnetic  dipole  aligned  with  the  Earth’s  rotation  axis  so 
that  the  magnetic  South  Pole  is  located  on  the  Earth’s  geographical  North  Pole. 

B  =  ^-\^\MT-Ur\xir-MT)  (3.1) 

471  r 

here  MT  is  the  Earth’s  magnetic  dipole.  The  strength  of  the  field  is 

B  =  _EL^_L^/1  +  3sin2^  (3.2) 

471  r 

where  X  is  the  latitude.  Using  polar  coordinates,  where  the  latitude  acts  as  the  polar 
angle,  we  obtain  the  polar  equation  for  the  force  lines 

r  =  R0  cos2  X  (3.3) 


Earth's  Dipole-like  Magnetic  Field 


Figure  0-1.  Near  Earth ’s  magnetic  field  topology,  according  to  the  dipolar 

approximation  considered. 
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R0  is  the  radial  distance  at  the  equator.  Using  this  relationship,  we  find  a  more 

convenient  expression  for  the  field  strength,  in  which  we  will  introduce  the  Mcllwain’s 
shell  parameter,  L,  which  at  the  magnetic  equator  corresponds  to  the  radial  distance  from 
the  Earth’s  center  expressed  in  units  of  Earth  radii,  and  labels  the  dipole  field  lines  (in  a 
multipole  field  that  is  no  longer  the  case). 

^  =  Po  Mr  -\/4-3cos2  X 
4n  cos6  X 

And  this  one  will  be  the  expression  used  within  our  calculation,  with 
MT  =  8.05 -1022  Am2. 

Background  Plasma 

We  will  consider  the  ionospheric  background  plasma  to  be  cold  and  with  only  two 
components,  electrons  and  protons. 

Density 

It  is  hardly  possible  to  give  an  analytical  expression  for  the  background  plasma  density  as 
a  function  of  its  dependences  on  the  various  parameters.  It  is  necessary  therefore  to  use 
empirical  values  or  numerical  models,  specifying  the  different  conditions  under  their  use. 
As  we  are  working  with  a  quasi-analytical  model,  we’ve  tried  to  model  the  plasma 
density  as  simple  as  possible.  Nevertheless,  the  modular  structure  of  the  code  makes  it 
easy  to  use  different  models  instead,  and  compare  the  different  results  obtained  with 
them. 


We’ve  taken  the  density  to  be  constant  along  field  lines,  which  is  consistent  with 
diffusive  equilibrium,  at  least  not  far  from  the  equator.  The  equatorial  density  has  been 
obtained  trough. 
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(3.5) 


where  z  is  the  altitude,  k  is  the  Boltzmann  constant,  Mo  is  the  ion  mass  (in  this  case  the 
proton  mass)  and  the  subscript  0  in  the  other  quantities  refers  to  their  value  at  a  reference 


altitude  of  zn  =  103  km . 
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Plasma  Density  in  the  Ionosphere 


Figure  0-2.  Plasma  density  as  a  function  of  the  distance  to  the  Earth ’s  center. 


Index  of  Refraction 


For  obtaining  the  index  of  refraction  of  the  background  plasma  we’ll  use  the  Appleton- 
Hartree  dispersion  relation,  that  is,  we  consider  the  plasma  to  be  cold,  neglecting  the  ion 
motions  but  considering  their  inertial  effects. 


2  B±F 

n  = - 

2  A 

The  definition  of  the  parameters  involved  is 

F 2  =B2 -4  AC 


(3.6) 


A  =  S' sin2  0  +  Pcos2  0 
B  =  RL sin2  0  +  Ps( 1  +  cos2  ©) 
C  =  PRL 
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And  they  are  functions  of  the  elements  of  the  dielectric  tensor,  which  for  a  2  component 
cold  plasma  are: 
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These  terms  depend  solely  on  the  wave  frequency  and  the  characteristic  frequencies  of 
the  plasma,  the  electron  and  proton  cyclotron  frequencies. 
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Emitter 


We  have  to  define  the  characteristics  of  the  waves  we  emit  in  frequency,  propagation 
angle  and  intensity.  To  do  so  we  have  first  to  determine  what  kind  of  emitter  we  will  use. 

Our  emitter  will  be  a  space  tether.  The  easiest  and  more  effective  way  of  using  a  space 
tether  as  an  emitter  is  making  it  work  as  a  half-wave  electric  dipole  antenna. 

The  orientation  of  such  an  antenna  must  be  perpendicular  to  the  local  magnetic  field  to 
maximize  the  emissions.  In  addition,  the  stiffness  of  the  tether  is  very  small,  making  it 
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necessary  to  look  for  a  stable  configuration  from  a  dynamical  point  of  view  to  avoid 
attitude  control  issues. 

We  have  chosen  as  configuration  that  of  a  tether  in  a  circular  equatorial  orbit,  an 
stabilized  by  gravity  gradient  (gravity  acts  as  a  restoring  force).  The  orientation  following 
the  vertical  in  an  equatorial  plane  keeps  the  tether  perpendicular  to  the  local  magnetic 
field,  and  the  circular  orbit  is  the  only  one  with  a  stable  equilibrium  position  for  the 
tether.  If  we  consider  the  Earth’s  magnetic  field  to  be  dipole-like,  the  antenna  will  work 
at  constant  L,  because  of  the  symmetry. 


Figure  0-3.  Tether’s  configuration  during  the  mission 


Now,  the  properties  of  an  electric  dipole  inside  a  plasma  will  define  the  wave  parameters 
needed  for  computation  of  the  diffusion  coefficient,  and  latter  the  characteristic 
precipitation  time. 

The  frequency  employed  is  a  design  parameter.  We  will  determine  it  among  the  permitted 
frequencies  trough  an  optimization  process,  being  the  optimization  criteria  the 
minimization  of  the  precipitation  time.  Once  we  define  the  frequency,  we  can  easily 
obtain  the  rest  of  the  parameters. 

The  first  one  is  the  tether  length.  The  antenna  efficiency  is  the  highest  when  its  length 
equals  half  the  wavelength  of  emission.  In  that  case,  the  efficiency  is  of  the  order  of 
«  70  -  80  % .  The  tether  length  is  thus  a  function  of  the  frequency,  the  propagation  angle 
and  the  background  plasma  conditions. 
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(3.7) 


2  nco 

The  propagation  angle  dependence  appears  in  the  refractive  index  of  the  background 
plasma. 


According  to  [Wang  &  Bell],  the  propagation  angle  value  for  a  perfect  electric  dipole 
inside  a  magnetized  plasma  is  only  a  function  of  the  frequency 
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and  this  will  be  the  value  we  use  in  our  calculations. 


(3.8) 


Finally,  the  emission  is  nothing  but  perfectly  monochromatic,  and  the  propagation  angle 
can  also  have  some  variation  around  the  value  given  by  (3.8).  We  will  use  typical  values 
for  their  dispersion,  the  same  we  can  find  in  the  literature. 

5(0  =  271-5  50  =  10°  (3.9) 

and  we  will  introduce  a  cut-off  at  low  frequencies 

®ic=K>m-5®  (3-10) 

Now  we  only  need  to  define  the  magnetic  field  intensity  of  the  wave,  Bw{k ).  If  we 

consider  an  electron  cyclotron  wave,  circularly  polarized  and  monochromatic,  which 
propagates  along  the  Earth’s  magnetic  field  lines,  we  can  apply  the  WKB  theory,  because 
all  the  magnetospheric  parameters  change  in  distances  very  large  in  comparison  with  the 
wavelength.  Under  this  assumptions,  the  magnetic  field  intensity  can  we  expressed  like 
this. 


i  cosfcot -  k{k)rdx\  +  j sinfoot - k(X)rdk 


(3.11) 


We  want  to  know  how  does  BW(X)  evolves.  On  one  hand,  it  varies  with  the  square  root  of 


the  refraction  index  of  the  background  plasma.  On  the  other  hand,  if  the  waves  are 
ducted,  which  is  another  assumption  we  are  doing,  the  intensity  varies  inversely  with  the 
surface  of  the  magnetic  duct.  Taking  into  account  both  effects,  we  have 


BM  =  K  o 


V  ^e,0n0 


(3.12) 
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The  factor  that  multiplies  the  intensity  of  the  magnetic  field  at  the  equator  is  of  order 
unity,  slightly  bigger  than  1  and  in  all  our  calculations  smaller  than  2.  We  can  then 
assume  the  intensity  of  the  wave’s  magnetic  field  to  be  constant,  and  equal  to  its 
equatorial  value.  Such  a  hypothesis  is  conservative,  because  it  assumes  a  smaller  value 
for  the  intensity,  which  yields  a  higher  precipitation  time. 

If  Bl  does  not  change,  it  is  just  a  multiplicative  constant  for  the  diffusion  coefficient  and 

the  diffusion  equation.  Since  this  square  is  proportional  to  the  antenna’s  emission  power, 
we  see  that  the  power  of  the  antenna  acts  solely  as  a  scale  factor.  Doubling  the  power  of 
the  system,  through  an  increase  in  the  power  of  the  satellites  or  an  increase  of  their 
number,  will  cut  the  remediation  time  in  half. 

To  find  an  expression  of  the  magnetic  field  intensity  as  a  function  of  the  emitted  power  is 
quite  delicate.  Following  the  scale  arguments  used  by  Inan  and  Bell,  we  have  taken  an 
intensity  of  15  pT  for  an  onboard  power  supply  of  10  kW.  We  will  use  this  value  for  our 
calculations,  but  we  should  not  forget  the  role  of  the  emitted  power  as  a  scale  factor. 

Numerical  Model 


Our  first  objective  is  to  obtain,  for  a  field  line  and  energy  given,  the  diffusion  coefficient 
as  a  function  of  the  particles’  equatorial  pitch  angle. 
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We  have  discretized  the  pitch  angle  interval  between  the  loss  cone  and  —  with  100 

points,  which  is  enough  for  all  the  cases  we  have  studied  at  a  reasonably  good 
computational  cost. 

To  discretize  the  problem  in  latitude,  we  have  taken  into  account  the  fact  that,  for  a 
certain  equatorial  pitch  angle,  the  weighting  functions  W;  are  zero  everywhere  but  in  a 
small  region  around  the  resonance  point.  Therefore,  we  have  calculated  first  for  each 
pitch  angle  the  resonance  latitude,  and  afterwards  we’ve  discretized  in  latitude  in  a  small 
interval  around  this  value.  The  amplitude  of  the  interval  depends  on  the  case.  After  each 
simulation  we  have  verified  that  the  amplitude  considered  was  enough  (zero  value  at  the 
boundaries,  unless  the  equator  is  one  of  them),  if  not,  it  was  widened.  The  number  of 
points  used  depends  on  the  dimension  of  the  latitude  interval,  ranging  from  100  to  200. 
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The  language  employed  to  code  the  program  is  Matlab®,  and  we  have  used  when 
possible  Matlab  functions  to  speed  up  the  calculations. 

Once  calculated  the  diffusion  coefficient  as  a  function  of  the  particles’  equatorial  pitch 
angle,  we  can  solve  the  diffusion  equation.  We  have  transformed  it  in  an  eigen  value 
problem,  with  matrix 


M  = 


0.3198-cosaf 


2D. 


1  + 


1 


sina 


o  J 


_ _  .  q  |  ^aa 

tan2a0  1.3802  -  0.31 98(sina0  +  -y/sina0 )  aa  da0 


da. 


+  D„ 


da. 


(3.13) 

We  have  discretized  this  matrix  using  finite  differences.  Before  choosing  the  type  of 
scheme  to  be  used,  we  need  to  have  a  look  at  the  order  of  magnitude  of  the  different 
terms  of  the  matrix.  If  the  second  derivative  is  the  most  important  term,  then  the  way  we 
discretize  the  first  derivatives  is  unimportant  from  the  numerical  stability  point  of  view. 
But  if  there  are  values  of  the  equatorial  pitch  angle  for  which  the  diffusion  coefficient 
tends  to  zero,  which  is  our  case,  the  first  derivatives  become  dominant,  and  we  can  no 
longer  neglect  their  influence  in  the  scheme  stability. 


The  stability  of  our  numerical  algorithm  depends  on  the  direction  taken  by  the 
derivatives.  To  decide  what  finite  differences  scheme  to  use,  we  have  done  a  boundary 
layer  analysis  to  evaluate  the  perturbations  evolution.  Through  this  analysis  we  have 
decided  to  employ  right  hand  finite  differences  for  the  first  derivatives.  For  the  second 
derivative  we  have  used  a  centered  scheme  of  second  order. 


The  problem  to  solve  is  then 


1  0  ••• 

.  o' 

gx 

gx 

0  •••  »v, 

mi  mM  •••  0 

gi 

-[*,]• 

gi 

0  . 

•••  0  -1  1 

_g n  _ 

_g n  _ 

(3.14) 


131 


where 
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T  —  ao,ic 

here  n  is  the  number  of  points  of  the  pitch  angle  discretization  and  A  =  — -  is  the 

n- 1 

integration  step.  This  scheme  does  not  give  any  numerical  problem  in  all  cases  studied. 
We  have  to  consider  also  =  Xn  =  0  because  of  the  boundary  conditions. 

Once  we  solve  this  problem,  the  eigen  values  give  the  characteristic  precipitation  time 
and  the  eigen  vectors  the  form  of  the  equatorial  pitch  angle  distribution  for  the  electrons, 

sK)- 
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Results 


In  this  chapter  we  will  analyze  several  cases,  so  as  to  decide  whether  or  not  such  a  system 
for  RBR  is  feasible.  First  of  all  we  will  decide  which  frequency  among  all  the  possible  is 
the  best  option  for  cleaning  the  radiation  belts  the  fastest  way.  Then,  we  will  analyze  the 
influence  of  the  distance  to  the  Earth  in  the  performance  of  the  system 

Frequency  optimization 

The  formalism  employed,  that  of  the  quasilinear  theory,  is  valid  in  a  certain  range  of 
frequencies.  The  upper  frequency  is  limited  by  the  absorption  (Landau  damping)  due  to 
the  background  plasma,  as  seen  in  §2.1.  We  have  taken  as  a  boundary  25%  of  the 
cyclotron  resonance  frequency.  The  lower  limit  corresponds  to  the  low  hybrid  resonance. 


to 


LH 


Q  Q. 


°y2  +  Q.A 

®„2+n«2 


(4.1) 


And  we  have  taken  as  a  boundary  twice  this  frequency.  Both  limits  define  a  range  in  the 
frequencies  we  can  use,  which  depends  on  the  distance  to  the  Earth. 


Limits  in  :he  frequencies  used 


Figure  0-4.  Possible  antenna  emitting  frequencies 
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To  analyze  the  influence  of  the  frequency  on  the  precipitation  time,  and  decide  which  one 
is  the  optimal  choice,  we  have  run  four  simulations  with  different  frequencies.  The  rest  of 
the  parameters  remain  the  same.  The  basic  conditions  for  all  of  them  are  L=1.5  and  an 
energy  of  1  MeV  per  electron. 

The  two  dominant  resonances  are  the  first  one  and  that  of  Landau.  We  have  plot  the 
diffusion  coefficients  corresponding  to  these  resonances,  together  with  the  latitude  of 
resonance  as  a  function  of  the  equatorial  pitch  angle.  The  four  frequencies  chosen  are  5 
kHz,  10  kHz,  15  kHz  and  25  kHz,  and  they  cover  the  whole  range  between  the  low 
hybrid  resonant  frequency  and  the  cyclotron  frequency. 


o 

g 

1  3 
b 
® 

< 

€  2 


-  5  kHz 
10  kHz 
15  kHz 

-  25  kHz 


Figure  0-5.  Left:  emission  frequency  influence  on  the  diffusion  coefficients.  Right : 
emission  frequency  influence  on  the  latitude  of  resonance 
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The  most  dramatic  changes  take  place  within  the  first  order  resonance.  The  more  we 
reduce  the  frequency,  less  particles  enter  in  resonance  with  the  wave,  but  the  higher  the 
diffusion  coefficient  value  in  at  the  loss  cone  boundary.  The  shape  of  the  diffusion 
coefficient  tends  to  become  more  flat  at  the  same  time.  If  we  plot  the  total  diffusion 
coefficient  as  a  function  of  the  equatorial  pitch  angle  for  the  four  frequencies,  we  obtain 

Diffusion  Coefficient: 


Figure  0-6.  Bounce  averaged  diffusion  coefficients  as  a  function  of  the  equatorial 
pitch  angle  for  different  values  of  the  emission  frequency.  In  all  cases  the  particles’ 

energy  is  1  MeV  and L=1.5 


Now  we  are  going  to  obtain  the  precipitation  time,  taking  into  account  all  resonances.  In 
fact,  we’ve  taken  into  account  just  the  first  5  resonances,  because  the  influence  of  the  rest 
is  negligible.  The  results  we  have  obtained  by  solving  the  eigen  value  problem  are  the 
following 
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Precipitation  time 
L=1.5 


Figure  0-7.  Characteristic  precipitation  time  for  the  electrons,  as  a  function  of  the 

emission  frequency 

The  characteristic  precipitation  time  has  a  strong  dependence  on  the  emission  frequency. 
We  will  then  try  to  use  frequencies  near  the  lower  limit  in  order  to  reduce  the  mission 
time. 

Distance  to  Earth  influence 

We  will  now  analyze  the  influence  of  the  distance  to  Earth  on  the  characteristic 
precipitation  time.  To  avoid  absorption  near  the  low  hybrid  resonance,  we  consider 
frequencies  equal  to  twice  the  low  hybrid  resonant  frequency.  We  have  calculated  this 
time  for  three  different  field  lines  and  for  a  frequency  near  twice  the  low  hybrid  resonant 
frequency. 


L  f  [kHz] 


1.5  10 
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2.0 


5 


2.5  2.5 

Table  4-1.  Frequencies  employed  for  different  distances  to  Earth 
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First,  we  will  see  the  effect  on  each  resonance  independently 
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Figure  0-8.  Left :  influence  of  the  distance  to  Earth  on  the  diffusion  coefficients,  for 
the  Landau,  first  and  second  order  resonances,  for  three  different  distances.  Right : 
influence  of  the  distance  to  Earth  on  the  resonance  latitude,  for  three  different 
distances  and  the  Landau,  first  and  second  order  resonances. 
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The  further  the  system  is  to  the  Earth,  the  resonances  take  place  nearer  the  mirror  points 
and  the  diffusion  acts  for  particles  until  higher  values  for  the  equatorial  pitch  angle.  As 
for  the  magnitude  of  the  diffusion  coefficients,  it  depends  on  the  resonance  we  consider. 


Diffusion  coefficient 


Figure  0-9.  Diffusion  Coefficients  for  three  different  distances,  a  frequency 
approximately  equal  to  twice  the  low  hybrid  frequency  and  particles  of  1  MeV. 


The  second  resonance  becomes  more  important  as  the  distance  to  Earth  increases.  The 
characteristic  times  obtained  when  solving  the  diffusion  equation  are 
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L  influence  in  Precipitation  Time 


Figure  0-10.  Characteristic  time  obtained  for  each  one  of  the  three  cases  studied 

previously 

There  is  no  evident  tendency,  but  all  of  them  are  of  the  same  order  of  magnitude,  and 
what  is  more  important,  they  are  short  enough  so  that  we  can  expect  this  to  a  be  a 
promising  RBR  method. 

Some  additional  antenna  considerations 

Taking  into  account  the  environment  conditions,  and  the  propagation  angle,  we  can  easily 
obtain  the  antenna  length  needed  for  each  situation  using  (3.7).  The  results  are 

L  Longitude  [m] 


1.5  164.3 

2.0  235.7 

2.5  293.4 


Table  4-2.  Antenna  length  needed  for  three  different  distances  to  Earth. 
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As  we  can  see,  the  lengths  needed  are  large  enough  to  make  unpractical  the  use  of 
conventional  dipole  antennas.  The  use  of  tethers  can  solve  this  problem.  However,  they 
are  short  enough  so  that  there  won’t  be  any  specific  production  or  control  problems.  The 
differences  in  length  are  significant,  thus,  if  we  wish  to  use  the  same  tether  to  clean 
different  L  shells,  we  will  have  to  design  it  so  as  to  have  a  variable  length,  to  avoid  severe 
efficiency  losses  when  we  deviate  from  the  initial  design  point.  The  adaptability  of  the 
system  depends  on  the  range  of  L  shell  each  tether  will  work,  which  depends  on  the 
number  of  tethers  and  the  area  to  be  cleaned.  There  is  then  a  trade  off  between  number  of 
satellites  and  complexity  of  each  one  of  them  and  the  design  will  be  chosen  in  terms  of 
cost  and  mission  time  (the  more  satellites  we  have,  the  less  time  it  takes  to  clean  a 
specific  area). 

Since  the  tether  length  is  small,  there  won’t  be  any  concerns  about  its  strength  and  for  the 
sake  of  simplicity  we  can  consider  a  design  with  constant  cross  section. 

The  tension  generated  through  gravity  gradient  is 


Figure  0-11.  Tether  configuration  and  nommenclature  for  the  calculation  ot  the  tether 

tension. 
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And  if  we  take  an  aluminum  tether  at  L=2  this  yields 


=  7.2  Pa 
5 


This  value  is  negligible  in  comparison  with  the  mechanical  strength  of  the  material.  We 
do  not  have,  then,  any  concern  about  the  system  construction. 

From  the  dynamical  point  of  view,  without  doing  a  deep  study  of  the  cable  dynamics,  its 
first  frequency  is  approximately  that  of  a  string  under  tension 

f\  =2A,rv®2'10”4^z  («> 

2 L  \  pi 

which  is  very  far  from  the  excitation  frequency,  which  is  that  of  the  waves,  because  the 
perturbation  forces  are  basically  due  to  the  interaction  between  the  Earth’s  magnetic  field 
and  the  currents  in  the  tether.  Thus,  there  should  not  be  any  attitude  control  problems  for 
the  system. 

Finally,  according  to  Edelman’s  research,  there  are  no  specific  sputtering  problems  due  to 
the  ions  bombardment. 


Synthesis 

If  we  take  as  a  characteristic  remediation  time  three  times  the  time  constant,  the  results 
that  we  have  obtained  are 


L  Remediation  time  [days]  Tether  length  [m] 


1.5  12  164 
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2.0 


9 


236 


2.5  15  293 

Table  4-3.  Results  for  the  different  distances  to  Earth 

In  this  analysis  we  have  just  considered  the  effect  of  YLF/ELF  waves.  There  are  other 
effects  that  affect  the  mean  life  of  a  trapped  particle.  This  other  effects  make  that  the 
mean  precipitation  time  for  a  1  MeV  particle  at  the  distances  we  are  considering  is  a  year 
[ Edelmann ].  Since  our  decay  times  are  small  in  comparison  with  those  obtained  without 
the  effect  of  VLF/ELF  waves,  we  can  infer  that  the  effect  of  the  latter  is  dominant,  and 
that  the  results  we  have  obtained  are  accurate. 

The  fact  that  the  decay  times  are  reduced  by  an  order  of  magnitude  or  more,  means  this  is 
a  promising  approach  to  Radiation  Belt  Remediation. 

We  have  not  taken  into  account  the  effect  of  possible  wave  reflections  in  the  ionosphere. 
Data  from  natural  whistlers  indicate  that  there  can  be  as  many  reflections  as  20  or  30.  We 
could  therefore  expect  further  reductions  in  the  remediation  time,  by  a  factor  of  10. 

To  take  into  account  the  effect  of  reflections,  we  would  need  an  in  depth  study  of 
VLF/ELF  waves  propagation  in  the  ionosphere,  which  is  out  of  the  scope  of  this  study. 


143 


Conclusion  and  recommendations 


Conclusion 

The  cleaning  of  a  thin  L  shell  by  a  10  kW  system  is  feasible  in  a  few  days  if  there  are  no 
particles’  sources  involved.  For  cleaning  significant  areas  of  the  ionosphere  we  need  a 
variable  length  system,  a  system  made  of  several  satellites,  or  a  combination  of  the  two, 
depending  of  the  mission  time  and  cost  we  are  targeting.  If  we  stay  relatively  near  the 
Earth,  the  RBR  time  stays  more  or  less  constant. 

This  method  appears  to  be  more  promising  than  the  one  studied  by  [Zeineh,  2005],  which 
consisted  in  a  high  voltage  tether  dispersing  the  ions  as  it  orbits  around  the  Earth.  The 
remediation  times  obtained  by  Zeineh  were  one  to  two  orders  of  magnitude  higher.  It  is 
true  that  for  comparing  both  methods  we  should  apply  the  Zeineh  method  to  electrons, 
but  this  needs  significant  modifications  of  the  prior  work 

The  results  obtained  depend  on  the  emitter  power.  If  we  multiply  the  power  of  the  system 
by  a  factor,  either  increasing  the  power  of  each  satellite,  or  increasing  the  number  of 
satellites,  the  total  remediation  time  will  be  reduced  by  the  same  factor. 

Finally  the  accuracy  of  the  results  can  be  improved  by  using  better  models  for  the  Earth’s 
magnetic  field  (IGRF,  ...),  the  ionospheric  plasma  density  (SLIM,  ...),  and  taking  into 
account  more  resonances.  The  modifications  needed  are  minor,  but  we  lose  many 
analytical  expressions. 

Recommendations 

Our  recommendations  for  future  research  are  the  following 

•  Development  of  a  code  for  calculating  the  ray  trajectory  of  YLF/ELF  waves  in  the 
ionosphere.  This  code  would  allow  us  to  better  know  the  area  of  influence  of  the 
antenna,  and  take  into  consideration  the  effect  of  multiple  reflections  near  the 
mirror  points.  Whistlers  can  reflect  several  times,  increasing  the  effect  over  the 
particles’  trajectories.  This  effect,  that  we  have  neglected  in  our  analysis  could 
reduce  even  more  the  remediation  time,  and  should  be  studied  in  detail.  Such  a 
code  could  be  based  on  the  WKB  optical  approximation  and  should  be  able  to  use 
several  different  ionospheric  models. 
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•  Numerical  simulation  of  the  VLF/ELF  emissions  by  an  electrodynamic  tether.  A 
PIC  code  specially  designed  for  the  study  of  electrodynamic  tethers  has  already 
been  developed  inside  SPL.  This  code  should  be  modified  to  allow  an  increase  in 
the  numerical  domain  without  increasing  the  computational  cost.  It  is  also 
necessary  to  introduce  boundary  conditions  that  allow  waves  to  cross  the 
boundaries  without  reflection.  Such  radiation  like  boundary  conditions  may  be 
similar  to  those  used  in  CFD  codes  for  aeroacoustic  simulations.  This  code  should 
allow  us  to  know  better  the  properties  of  the  waves  generated  by  the  tether, 
improving  then  our  results. 

•  Numerical  study  of  the  wave-particle  interaction.  The  region  in  which  the 
resonance  takes  place  is  small  in  comparison  with  the  particles’  trajectories  (under 
6  degrees  in  latitude)  but  too  large  to  simulate  the  resonance  phenomenon  through 
a  PIC  code.  If  we  consider  that  the  wave  is  not  perturbed  by  the  particles,  we  can 
assume  they  are  not  coupled  and  we  can  try  a  “test  particle”  approach  to  the 
problem.  Then  using  statistical  Monte  Carlo  techniques,  we  can  obtain  the 
average  pitch  angle  scattering  and  the  precipitating  particle  fluxes,  comparing  the 
results  with  those  obtained  with  the  quasilinear  theory  with  a  reasonably 
computational  cost. 

•  Study  of  the  influence  of  ions  other  than  protons  inside  the  ionosphere  and  of 
other  kinds  of  wave  (EMIC,  Chorus). 

All  these  studies  and  simulations  should  allow  us  to  know  better  the  behaviour  of  such  a 
system  before  starting  with  flight  models. 
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Emission  of  Whistler  Waves  from  an 
ionospheric  tether 


May  26,  2008 


1  Equations 


An  electromagnetic  wave  in  plasma  is  formulated  by  application  of  Maxwell’s 
equations  and  the  equations  of  motion  of  the  particles.  In  our  problem,  we 
use  Maxwell’s  equations  and  the  electron  equation  of  motion.  We  neglect 
the  ion  motion  because  the  frequencies  of  interest  are  above  the  lower  hybrid 
frequency. 

Maxwell’s  equations  (Neglecting  displacement  current): 


V  x  = 


dB{ 


(C) 


dt 


V  x  B(a  = 


(r)  —  ftoJ(?) 


Electron  equation  of  motion 


(1) 

(2) 


me 


dve 

dt 


e(£)v)  +  ve  x  i?0(c))  —  meueve  (  ue  :  collision  frequency) 


(3) 


From  Eq.(l)  andEq.(2), 


2  dve 

V  -E(C)  V( V  *  E(C))  —  ft (W^eO  \~  ft o 


djs 

dt 


(4) 


147 


2  Fourier  Transformation 


We  can  solve  those  equations  using  Fourier  transformation.  In  this  problem, 
we  analyse  the  wave  propagation  with  fixed  frequency  oj. 

One  of  the  space-time  Fourier  components  can  be  expressed  as 


E(rl  =  E(u,k)  exp  l(LOt-k-x) 


From  Eq.(4), 


—k2E  +  k(k  ■  E)  —  —i/JozneoU?ve  +  ijio^js 


From  Eq.(3), 

me(iuj  +  ue)ve  —  —e(E  +  ve  x  B0) 
Then,  from  Eq.(6)  and  Eq.(7), 


E 


1  - 


ik2  +  ve) 

jjL0u  e2ne0 


ik 2 


/i0ne0ue 


E  x  Bq  +  (k  •  E)i 


(5) 

(6) 
(7) 

k  me(iuj  +  ve) 
jjL0u  e2ne0 


-\ 


-me(iu  +  ue)  7  js  x  B0 

~Js 


e2ne0 


cneQ 


Dehne 


bo  —  Bo/ Bo 

s=*- 


v 


LOt 

CO 

OJ, 


ce 


ce 


1  = 


mP 


c  neofio 
I<  =  kl 


®ce 


e2ne0 


(9) 

(10) 

(11) 

(12) 

(13) 

(14) 


Assuming  that  B0  =  10-5  Tesla  and  electron  temperature  is  1  eV,  the  plasma 
skin  depth  l^ace  and  8  for  n  =  lO10  [/m3]  and  n  =  108  [/m3]  are  estimated  as 
shown  in  Table  2.  (toce  ~  1.7  x  106  [rad/sec]) 


k  x  B0 
li0ene0  to 

(8) 
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Table  1:  Estimated  l7crce  and  8 


n  [/m3] 

l  [m] 

ace  [Si/m] 

5 

1040 

50 

1.6  x  10"4 

4  x  10"8 

108 

500 

1.6  x  10"6 

4  x  lO"10 

Then,Eq.(8)  becomes 
E  ( 1  +  K 


2  J<28\  .K2  £  A/  ^  -^x6o 

—  % -  —  % — E  x  60  T  (A  •  E)  —K  +  % - b  i—K 

V  j  V  \  /  \  v  V 


(iu  +  J)  js  +  js  X  bo 


(15) 


The  configuration  is  such  that  k  is  along  x-axis, and  B0  makes  an  angle  9  to 
k  and  is  in  the  x-y  plane  (Fig.(l)). 


Z 


y 


Figure  1:  k  and  B0  configuration 
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Then, 


bo  —  cos  Ox  +  sin  By 

E  x  bo  —  —  sin  0Ezx  +  cos  0Ezy  +  (sin  0EX  —  cos  OE^j  z 
K  x  bo  —  K  sin  Oz 

js  -  -  sin  Ojszx  +  cos  Ojszy  +  (sin  0jsx  -  cos  Ojs^j 
(£,  y,  unit  vectors  on  x,  y,  z  axes) 

Thus,  Eq.(15)  becomes  the  system 


-  .A2  .  -  (iu  +  5)  jsx  —  sin9jsz  ,  , 

Ex  +  i — sin  9EZ  =  -1 —  u  - —  (16) 

V  (Joe 

C  K2  .K2S\  p  .K2  p  (iu  +  6)  jly- cos  0fsz 
1  +  A  —  i -  Ev  —  i cos  0EZ  — - 17 ) 

V  v  )  v  (Tce 

■  K2  ,  (,  ,  Ts2  -K28\  p,  (iu  +  8)jsz  +  sm0jsx-cos0jsy 

% - cos  vEy  +  1  +  A  —  % -  Ez  = - 

v  \  V  (JCR 


Ey  and  Ez  are  decoupled  from  Ex  (along  fc),  but  Ex  is  tied  to  Ez  by  Eq.(16). 


3  Dispersion  relation 

The  determinant  of  Eq. (16), (17), (18)  is  same  as  the  determinant  of  Eq. (17), (18). 
For  a  valid  EXJ  Eyj  EZJ  we  need  this  determinant  to  be  zero. 

('1+A-2_I£hy_iico^=0  (19) 

\  V  J  V1 

And  the  roots  are 


These  are  the  dispersion  relations.  As  in  Fig. (2),  for  9  <  cos  1  v,  Ku  repre¬ 
sents  slightly  damped  waves,  and  Ah  represents  heavily  damped  waves.  For 
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cos-1  v  <  9  <  7T  —  cos-1  v,  Ku  and  Kd  represent  heavily  damped  waves.  For 
9  >  7T  —  cos-1  i/,  Ku  represents  heavily  damped  waves,  and  Kd  represents 
slightly  damped  waves.  So,  the  Ku  wave  propagates  when  9  <  cos-1  and 
the  Kd  wave  propagates  when  9  >  tt  —  cos-1  v. 


Figure  2:  Ku  and  Kd  pole  position 


4  Matrix  form  of  Ohm’s  law 

Eq. (16), (17), (18)  can  be  solved  as  follows 


Ex 

3  sw 

E  = 

Ey 

N 

II 

,  i*s 

N 

II 

3  sy 

.  Ez  _ 

_  jsz  _ 

(22) 


Z 


iu  (1  +  K 2)2  —  (A"2  +  sin2  9) 

i—  sin  9  cos  9 
(l  +  AT2)sin.0 


sin  9  cos  9 

iu  (1  +  K 2)  —  (cos2  9) 
—  cos  9 


—  (1  +  K 2)  sin# 
cos  9 

ip  (1  +  K 2)  —  (cos2  9) 


(23) 
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2 


(24) 


A  = 


1  +  K2 


,K2S 


V 


5  Antenna  related  axes 


For  the  computation  of  the  inverse  Fourier  transformation,  we  will  need  to 
vary  k  in  all  possible  ways.  Thus,  it  is  better  to  change  the  coordinates  to  ones 
which  do  not  change  in  the  inverse  Fourier  transformation  process.  Now,  k  is 
on  the  x-axis.  We  will  change  the  coordinates  by  defining  axes(Wjg,  Yg,  Zb) 
such  that  Bo  is  on  the  z-axis  (B-frame)  as  in  Fig.  (3).  The  antenna  is  in  the 
Xb~  Zb  plane  and  has  an  angle  a  to  Bo.  k  is  oriented  according  to  the  polar 
angles  (#,</>). 

Any  vector  can  be  expressed  in  the  old  coordinates  V  =  (Vx^Vy^Vz)  and  in 
the  new  coordinates  Vb  —  (Wb,  Fy#,  Vzb)  and  the  two  representations  are 
related  as 

Vb  =  RV  (25) 

R  is  a  rotation  matrix  and  has  components  as  follows 


R  = 


sin  9  cos  </>  —  cos  9  cos  </>  sin  </> 
sin  9  sin  </>  —  cos  9  sin  </>  —  cos  </> 

cos  9  sin  9  0 


(26) 


R-1  =  RT,  so  V  =  RtVb. 

From  Eq.(22),  the  current  to  held  relationship  is  E  —  Zjs.  Transforming  to 
the  B-frame  gives 


and,  therefore 


RZj5.  =  RE  —  Eb  —  ZbJsB  —  ZbRJs 
Zb  =  RZRt 


(27) 

(28) 
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Zb 


Figure  3:  Antenna  Related  Axes 


Explicitly,  Zb  Fas  components  as  follows 


i 

<J  ceA 


— p  +  K2  sin2  9  (—p  cos2  (j)  +  ^  sin2  (fj 
—  (l  +  K2  sin2  9 )  —  iK2  (1  +  K2)  (^  —  v)  sin2  9  sin  <j>  cos  <f> 
—K2  ( p  cos  (j>  +  sin  </>)  sin  0  cos  0 

—  (l  +  K2  sin2  9 )  —  i  AT2  (1  +  K2)  p  —  i/)  sin2  #  sin  <f>  cos  <f> 
—  (p  —  sin2  cos2  <j>  —  p  (l  +  A"2  sin2  0)  sin2  <^> 

A"2  (—  p  sin  o  +  cos  o)  sin  0  cos  0 


A"2  sin  0  cos  9  (sin  <f>  —  p  cos  o) 
—A"2  sin  9  cos  9  { p  sin  cj>  +  cos  <^>) 
— pK 2  cos2  #  +  iv  (1  +  K 2) 


(29) 


where 


—  iz/  (l  +  A"2) 


(30) 
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6  The  antenna  model 


For  an  antenna  source  current  distribution  js{x),  the  Fourier  transformation 
is 


1 

(2^7 


js(x)el3xd3x 


(31) 


and, 


Js{x) 


a  inside  wire 

Trar 

0  outside  wire 


(32) 


a  is  the  radius  of  the  antenna,  la  is  the  unit  vector  along  the  antenna,  Is  is 
the  current,  za  is  distance  along  the  antenna.  So,  if  /0  is  Is  (za  =  0),  Eq.(31) 
becomes 


1  To  f 

(2  nfna2 


h  (za)  ik-g  1 3 


<Tx 


antenna 


(33) 


Let  k_ l  be  the  projection  of  k  on  the  plane  perpendicular  to  the  antenna, 


Figure  4:  Antenna  and  k 

and  k\\  be  the  projection  of  k  on  the  antenna.  Then  if  kj_  makes  an  angle  x  to 
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the  xa  axis  (which  is  on  the  XbYb  plane),  and  x  has  cylindrical  components 
(R^,za)  in  (xa,ya,za), 

k  •  x  —  k^Rcos  (ip  —  x)  +  k\\za  (34) 

The  integration  volume  element  is  d3x  =  RdRd^dza,  Using  Bessel  functions 
(dn(x)  =  ^  elxsin0~in0  d9^j ,  the  integral  on  ^  is 


* 2 7T 


eik±R***-x)d^  =  27rJ„(k±R) 


(35) 


and  the  integral  on  R  is, 


2  7T 


/  2ttJo  (fcj_i2)  7U1R  =  -j  (akj_)  Ji(ak±) 
I  o  k±_ 


(36) 


Then,  we  can  get 

h  (k)  = 


1  I0  2na  T  /  7  N  fL/2  Is{za)  %k 


(2  ,fXal^Jl{akL) 


-L/2  h 


e8fcH  Zadzn 


(37) 


Integrate  by  parts  and  assume  Is(za  =  ±L/2)  =  0,  and  U(^a)  =  Is(-za) 


^  IS(Zg) 

-L/2  /0 


2  rL'2  d  fi 

e  = -fa) /,  sin(^W  UW  (38) 


31  ^0 


We  next  assume  a  triangular  current  distribution  for  Is(za ): 


C  _  i  _  2j£a 

/o  “  L 


Then, 


r'L'2  h(Za) 


.1/2  h 


Substituting  this  in  Eq.(37), 


(39) 


(40) 


(41) 
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(42) 


In  general,  ak±_  <  1,  so  Ji(ak_\_ )  ~  |afc_ |_,  and  then 


4I0L  1  -  cos  (hL/2) 
h[k)~(2  tt)3  k\I? 


For  long  waves  (^1—  -C  l) 


jAk)  =  a 


(43) 


2  (27t)3 

This  is  the  Fourier  transformed  delta  function  at  the  origin  with  strength 
2-fo-i 


UoF. 


7  Inverse  Fourier  transformation 


We  now  have  explicit  expressions  for  Z  and  js ,  and  the  Fourier  transformed 
E  held  is  E  —  Zsjs-  Then,  the  inverse  Fourier  transformation  is 


E  (x) 


Elk 


—  ik •; 


:d3k 


(44) 


E  —  ZsAjs  with  A  — 


sin  a 

0 

cos  a 


Zb  is  shown  in  Eq.(29) 


(45) 


d3k  =  k2dk  sin  9d9d<$> 


The  process  for  integration  is  as  follows. 

1.  Fix  #,</>,  and  integrate  on  k  from  0  to  oo.  Integration  on  k  is  effected 
in  the  complex  k  plane. 

2.  For  9  and  </>  integration,  we  will  use  the  stationary  point  method  (ex¬ 
planation  in  [10  Stationary  point  method]).  After  k  integration,  for 
large  kr  (far  held)  identify  the  values  (0,  <f>)  that  make  the  phase  (k  •  x) 
stationary  for  each  observation  point  x. 
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3.  We  can  carry  out  the  9 ,  </>  integration  by  covering  the  vicinity  of  the 
stationary  points  only.  At  the  stationary  point,  (0S,</>S),  we  can  take 

the  E  term  out  of  the  integral  and  also  we  can  expand  (k  •  x)  to 
quadratic  order  in  9  —  9S  and  </>  —  <f>s. 

4.  Perform  the  integration  of  e~lk'x  (k  •  x  has  only  zeroth  and  second  order 
terms  now).  This  completes  the  process  for  most  observation  points  x. 

5.  For  certain  observation  points  x  where  the  second  derivative  of  k  •  x  is 
zero,  the  expansion  must  be  carried  out  to  the  next  order. 

For  observation  points  along  B7  we  need  another  treatment  for  the 
integration  because  the  <f>  derivative  of  k  •  x  is  zero  at  that  point. 


Figure  5:  Antenna,  observation  point  and  k 
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8  Causality  considerations 


At  some  distance  from  the  antenna,  information  and  power  travel  at  the 
group  velocity. 

*  =  |  W 

We  define  polar  angles 
#,  </>  :polar  angle  of  k 


9X7  4>x  :polar  angles  of  the  vector  to  the  observation  point 
9q7  (j>o  :polar  angles  of  vq 

The  direction  of  vq  should  be  same  as  the  direction  of  the  observation  point. 
So,  for  the  waves  to  be  observed  from  (9X,  (j>x), 


9cj  —  <f>G  —  4>x 


(47) 


should  be  satisfied.  For  very  small  damping  (d  -C  //),  the  dispersion  relations 
(Eq.(20)  and  Eq.(21))become 


So,  we  define 


K2  = - - - 

“  cos  0  —  V 


or  u  — 


K2  cos  9 

i  +  K 


-v 


Kd  =  -  n 

cos  9  +  u 


or  v  — 


K\  cos  9 

'TTW 


&u- 

&d+ 

Kel¬ 


ts 


cos  9  —  u 

i  * 

v  cos  9  —  u 


—u 


cos  9  +  u 


—u 


cos  9  +  u 


(48) 

(49) 


(50) 
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A„  propagates  only  for  0  <  cos-1  //,  Kd  propagates  only  for  ()  >  ~  —  cos-1  //, 
where  //  =  -j-.  Using  Kx  —  A  sin  d  cos  o,  A ;/  =  A"  sin  #  sin  </>,  Ah  =  A  cos  A 
so  for  A  „  wave, 


Vqxu  — 

^Gyu  — 

VfJzu  — 


VGx 

zJ  ,:r  / 


cA  1 
0I\-x 
duj  1 


*9  Ay 


VGz 


du  1 


dKz  u}c^l  1  A  A 


1  -  K2  Kx 
~  Xz(l  +  K2f-K 
..  1  -  A2  A, 

^  (1  +  A"2)2  A 
K  f  1  -  K2 
+  1  +  A"2 


cos2  9 


and  for  the  Kd  wave, 


Vcjxd  —  —A 
Vbyd  =  —  A 
Vgw  = 


1  -  A2  Kx 
(1  +  A"2)2  "A 
1  —  A"2  Ah 


(1  +  K2)2  K 


K 


1  +  K2 

From  Eq.(50),  Eq.(51),  and  Eq.(52), 


1  + 


1  -  K 
1  +  K 


cos2  9 


COS  Au+  = 


COS  9qu—  — 


cos9Gd+  - 


cos  9Gd-  - 


2  cos  0  (cos  9  —  u)  +  sin2  0 
-y/ 4  (cos  9  —  u)2  +  sin2  9 
2  cos  9  (cos  9  —  u)  +  sin2  9 
\J 4  (cos  0  —  u)2  +  sin2  0 
2  cos  0  (cos  0  +  v)  +  sin2  0 
\j 4  (cos  0  +  i/)2  +  sin2  0 
2  cos  9  (cos  9  +  v)  +  sin2  0 
-v/ 4  (cos  0  +  i/)2  +  sin2  9 


(51) 


(52) 


(53) 


There  are  three  methods  to  visualize  the  wave  propagation: 

1.  VG  graph 

2.  k  tip  graph 
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3.  9  vs  9%  graph 

The  Vcj  graph  is  a  locus  of  the  components  of  Vcj  along  B0  and  perpendicular 
to  Bo-,  with  the  k  direction  9  as  a  running  parameter.  Fig. (6)  is  the  Vcj  graph 
for  //  =  0.1.  For  the  Ku  wave,  Vcj  is  on  the  5  axis  for  9  —  0,  then  9cj  increases 
until  it  reaches  a  maximum,  then  it  decreases  to  zero  at  9  —  cos-1  (2 v)  and 
goes  to  negative.  When  9  —  cos-1  u7  the  Ku  wave  is  resonant.  The  Ku  wave 
doesn’t  propagate  for  9  >  cos-1  v,  so  that  branch  ends  at  this  point. 


Figure  6:  Vcj  graph 
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The  k  tip  graph  is  a  locus  of  k  with  9  as  a  running  parameter.  Vcj  is 
perpendicular  to  the  locus  of  k  at  each  point.  Fig. (7)  is  the  k  tip  graph  for 
v  —  0.1.  For  the  Ku  wave,  9  cannot  be  larger  than  cos-1  so  the  k  tip  curve 
for  Ku  goes  to  the  asymptotes  9  —  cos-1  v. 


o 


180 


Figure  7:  k  tip  graph 
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Using  Eq.(53),  we  can  draw  the  9  vs  9X  graph  because  for  the  wave  to 
propagate,  9X  should  be  equal  to  9q.  Fig. (8)  is  the  9  vs  9X  graph  for  v  —  0.1. 
For  one  observation  point  with  9X,  only  the  waves  with  9  on  those  curves  in 
the  graph  can  propagate.  In  Fig. (7),  Fig. (8)  and  Fig. (8),  @,  (b),  (c)  and  @ 
corresponds  to 

(a):9  =  0,  @:9q  max,  (c ):9  =  cos-1  2 @:9  =  cos-1  u 
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The  branch  that  is  observed,  as  well  as  the  particular  k  within  that 


branch,  changes  as  the  observation  point  angle  9X  changes.  First,  if  the 
observation  point  is  very  near  to  the  5  axis,  which  means  9X  &  0  but  pos¬ 
itive,  three  waves  (2  Ku+  waves  and  1  Kd-  wave)  propagate  as  in  Fig.  (9), 
Fig. (10),  and  Fig. (11)  at  (f>  —  cj>x.  Also  three  waves  (2  Kd-  waves  and  1  Ku+ 
wave)  propagate  to  the  points  very  near  to  the  5  axis  at  </>  =  tt  +  only 
”  behind”  B$. 


VGz 


0.6 


0*4 


0.3 


0.2 


0.1 


-0.6 


-0.5 


-0.4 


-0.3 


-0.2 


-0.1 


-0.1 


-0.2 


3  waves 


-0.3 


-0.4 


-0.6 


Figure  9:  Region  1  Vg  graph  (Propagation  along  Bo) 
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0 


Figure  10:  Region  1  k  tip  graph  (Propagation  along  Bq) 
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Figure  11:  Region  1  9  vs  0X  graph  (Propagation  along  Bo) 
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In  region  2,  where  0  <  9X  <  sin-1  three  waves  (2  Ku+  waves  and  1  K(j- 
wave)  propagates  as  in  Fig. (12),  Fig. (13),  and  Fig. (14).  Again,  there  are 
three  other  waves  observable  from  the  (nearby)  symmetric  direction  behind 

Bo. 


Figure  12:  Region  2  Vq  graph  (very  near  B0) 
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0 


Figure  13:  Region  2  k  tip  graph  (very  near  Bq) 
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Figure  14:  Region  2  9  vs  0X  graph  (very  near  Bq) 
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On  the  point  where  9X  —  sin-1  three  waves  (2  Ku+  waves  and  1  A,/_ 
wave)  propagates,  of  which  the  K(i~  wave  is  resonant  as  in  Fig. (15),  Fig. (16), 
and  Fig. (17).  From  Fig. (16),  we  can  easily  know  that  9X  at  the  resonant 
point  is  9X  —  |  —  cos-1  v  —  sin-1  v 


Figure  15:  Region  3  Vcj  graph  (on  the  resonance  cone) 
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Figure  16:  Region  3  k  tip  graph  (on  the  resonance  cone) 
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Figure  17:  Region  3  0  vs  (),.  graph  (on  the  resonance  cone) 
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In  region  4,  where  sin  1  v  <  9X  <  cos  1  —  ^  u  '  two  waves  (2 

Ku+  waves)  propagates  as  in  Fig. (18),  Fig. (19),  and  Fig. (20). 


Figure  18:  Region  4  Vcj  graph  (outside  resonance  cone,  but  within  propaga¬ 
tion  cone) 


172 


0 


Figure  19:  Region  4  k  tip  graph  (outside  resonance  cone,  but  within  propa¬ 
gation  cone) 
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Figure  20:  Region  4  9  vs  9X  graph  (outside  resonance  cone,  but  within  prop¬ 
agation  cone) 
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On  the  cone  where  9X  —  9X max  =  cos" 


4(l-^2) 


,  only  one 


wave  (1  Ku+  wave)  propagates  as  in  Fig. (21),  Fig. (22),  and  Fig. (23).  9q 
reaches  its  maximum  and  no  wave  can  propagate  in  the  region  where  9X  > 


cos 


-l 


3V2(  l-^)-2uy/^ 


.  From  Fig. (23),  |^f-  =  0  on  this  point.  This  yields 


From  Eq.(53), 


cos  9  —  v  + 


1  -  u2 
3 


9q  —  cos 


3 


4(1  -v2) 


(54) 


(55) 
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VGz 


Figure  21:  Region  5  Vg  graph  (on  the  edge  of  the  propagation  cone) 
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Figure  22:  Region  5  k  tip  graph  (on  the  edge  of  the  propagation  cone) 
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Figure  23:  Region  5  9  vs  0X  graph  (on  the  edge  of  the  propagation  cone) 
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The  same  (symmetric)  discussion  applies  in  the  vicinity  of  9X  —  tt,  this 
time  involving  the  Ku-  and  Kd+  branches. 


9  K  integration 


As  a  first  step  in  computing  the  inverse  Fourier  transform  (Eq.(44)),  we  fix  the 
observation  point  x(r,9x,<f>x)  and  the  wave  direction  (#,</>)  of  k  —  K/l,  then 
perform  k  integration.  Now,  allow  K  to  be  a  complex  variable  and  perform 
complex  integration  using  the  residue  theorem,  although  the  K  integration  is 
in  principle  along  the  real  axis.  In  order  to  use  the  residue  theorem,  we  need 

to  close  the  path  of  integration  around  some  of  the  poles  of  A  (A"),  which  are 
Ku±  and  Kd±  as  given  by  Eq.(20)  and  Eq.(21).  This  closure  should  be  done 
in  the  half  plane  in  which  the  e%k'x  does  not  diverge.  Thus  the  imaginary 
part  of  k  •  x  should  be  negative  in  this  closure. 

Now 


k  •  x 


T 

—  K  [sin  9  sin  9X  cos  (</>  —  <f>x)  +  cos  9  cos  9X\ 

r  Tr 

-A  cos  7 


(56) 

(57) 


7  is  the  angle  between  k  and  T,  and  cos 7  is  the  bracket  in  Eq.(56).  Thus, 
when  cos  7  >  0,  we  must  have  lm(K)  <  0,  and  vice  versa.  Since  vq  should 
be  along  T,  7  is  also  the  angle  between  k  and  vq. 

We  limit  discussion  to  Ku+  and  Kd- ,  because  the  Ku-  and  Kd+  case  follows 
by  symmetry.  From  the  k  tip  graphs  in  the  causality  considerations  section, 
the  angle  between  k  and  vq  is  always  less  than  90°.  So,  cosy  is  always 
positive.  Then,  the  integration  path  and  relevant  poles  are  shown  circled 
in  Fig.  (24).  With  the  residue  theorem,  the  inverse  Fourier  transform  is  now 
reduced  to  a  double  integral  as  follows. 


E{f)  ~  0 


27ibRes 


4> 


k2 E^etk'x@poleKU7d  sin 9d9d<$>  (58) 


(the  |  is  because  the  original  K  integration  is  in  (0,  00)  ,  not  (  —  00,00)). 
Poles  are  from  T  'm  Eq.(29)(AA  and  Kd  are  the  roots  of  A  =  0).  Thus,  the 
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Figure  24:  K  and  integration  path 


residue  of  ^  is 


Res 


A 


K  -  K 


u± 


COS^ 


(K2  -  IQ)  (IO  -  Kl) 


=F (cos  9  —  v) 


A  — A 


4  cos  9 


for  K,,± ,  and 

1 


Res 


A 


K  -  K 


d± 


COS^ 


(K>-Kj)(K>-Kl) 


(59) 

± yV  (—  cos  9  —  i/) 


K-Kd± 


4  cos  9 


(60) 

for  Kd±. 

We  can  get  Res  E^e‘k'x@pol eKU}d  using  Eq.(59),  Eq.(60),  and  substituting 

A  „  and  Kd  into  Eq.(44)  and  Eq.(45).  Then,  considering  the  direction  of  the 
integration  paths,  the  integration  is 


E(g)  — 


f  f  2iri  (  Zb  j  Aajse%  sin  9d9d(f>  (61) 

'*  Vvalid  waves  / 
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sin  a 


(62) 


ZB{KU± ) 

— p  +  K 2  sin2  0  (— p  cos2  <j>  +  ^  sin2  </>) 

—  (l  +  A"2  sin2  0)  —  i K2  (1  +  K2)  (4  —  u)  sin2  9  sin  cf>  cos  cf> 
—K2  ( p  cos  cj>  +  sin  <f)  sin  9  cos  9 

—  (l  +  K2  sin2  9 )  —  iK2  (1  +  K2)  (4  —  u)  sin2  9  sin  <f>  cos  <f> 

—  (yP  —  sin2  9^j  cos2  <f>  —  p  (l  +  AT 2  sin2  #)  sin2  <f> 

K 2  (— p  sin  cj>  +  cos  </>)  sin  0  cos  0 
A"2  sin  0  cos  9  (sin  cj>  —  p  cos  <f) 

—K2  sin  9  cos  9  (p  sin  cj>  +  cos  <f)  (64) 

—pK2  cos2  9  +  iu  (1  +  K2)  K_K 


ZB(I<d±) 

—p  +  K2  sin2  9  (— p  cos2  </>  +  ^  sin2  </>) 

—  (l  +  AT 2  sin2  0)  —  i AT2  (1  +  A"2)  (4  —  z/)  sin2  0  sin  <j>  cos  <j> 

— K 2  (p  cos  cj>  +  sin  </>)  sin  0  cos  0 

—  (l  +  A"2  sin2  0)  —  i AT2  (1  +  K2)  (4  —  z/)  sin2  0  sin  <f>  cos  <f> 

—  (p  —  sin2  9^j  cos2  <j>  —  p  (l  +  AT 2  sin2  0)  sin2  <j> 

K 2  (— p  sin  cj>  +  cos  <f)  sin  9  cos  0 
K2  sin  9  cos  0  (sin  cj>  —  p  cos  <f) 

— K2  sin#  cos  9  (psin  cj>  +  cos  <f)  (65) 

—pK2  cos2  9  +  iu  (1  +  A"2) 

v  7  J  K—Kd+ 


10  The  stationary  point  method 

For  the  9  and  </>  integrations,  we  will  use  the  stationary  point  method.  For 
the  integration  Eq.(61),  if  r  is  large,  which  means  the  observation  point  is  far 
from  the  antenna,  the  exponential  term  oscillates  much  more 
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than  the  other  terms  \  27riZ^u  ^  Ajs  sin  9j  .  However,  at  some  points  where 

If  =  If  =  05:  the  exponential  term  does  not  vary  rapidly.  These  points  are 
called  stationary  points.  For  example,  Fig. (25)  is  the  real  part  of  the  inte¬ 
grand  (x  component)  with  r  =  1000km.  There  are  two  stationary  points  in 
this  graph  (red  dots).  For  wave  orientation  (#,</>)  not  near  these  stationary 


Figure  25:  Integrand 

points,  the  integrand  (E  e~ik'x)  contribution  to  the  integration  is  small, 

because  e~lk'x  term  oscillates  rapidly,  and  cancellations  occur  in  the  integra¬ 
tion. 

So,  near  the  stationary  point  (9SJ  </>5),  we  can  take  ^2t \iZ(Ku  d)Ajs  sin terms 
out  of  the  integral  and  we  only  have  to  integrate  the  exponential  term.  In 
addition,  since  only  the  vicinity  of  (0SJ4>$)  contributes,  the  exponent  can 
be  expanded  in  a  Taylor  series  about  this  point,  where  the  first  derivatives 
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are  zero,  and  terms  higher  than  second  order  in  (9  —  9S)}  (</>  —  cf>s)  can  be 
neglected.  Fig. (26)  is  the  plot  of  the  stationary  points  =  ||-  =  0^  for 

observation  direction  9X}<f)x.  The  propagating  wave  lines  in  Fig. (8)  are  ex¬ 
actly  on  these  stationary  points  lines.  So,  in  the  stationary  phase  method,  the 
waves  which  do  not  propagate  to  the  observation  point  are  cancelled  mathe¬ 
matically,  and  only  the  waves  which  propagate  to  the  observation  point  (they 
are  stationary  points)  contribute  to  the  integration. 


Figure  26:  Stationary  points 
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11  Propagation  near  B 

If  the  observation  point  is  very  near  the  magnetic  held  line  from  the  antenna, 
there  is  a  interesting  phenomenon  about  wave  proapagation.  For  this  analy¬ 
sis,  9X  in  Fig. (27)  needs  to  be  very  small.  After  K  integration,  Eq.(61)  can 


Figure  27:  Observation  point  near  B  line 


Wtt 


be  written  as 

E  =  %7viG(es^x)  /  e  sinOdO  /  e  ^ 

^  Jo  Jo 

with 


(J(QS$X)  — 

And,  using  a  Bessel  function, 


E 


,4  =  7'/" 


E  -  47t  Jo  [A 


sin  9  sin  9r 


X  \  —%A 


a /  cos  9  —  v 


e  Vcose-F  sin 


(66) 


(67) 
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If  the  Jo  function  can  be  moved  out  of  the  integral  (this  is  valid  only  when 
^<1,  see  discussion  below)  we  obtain 


E 


4ff2iG(es,<f>x)Jo 

4ff2iG(es,<f>x)Jo 


sin  9S  sin  9X 
A — , 

V  cos  9  —  u 

sin  9S  sin  9X 
A — , 

V  cos  9  —  v 


cog  9  cog  &x 
g  y/cos  9  —  v 


sin  9d9 


e~iA®'  sin  9d9  (68) 


^ ,  cos  9  cos  9X 
= 

V  cos  9  —  u 

The  stationary  points  for  are  9  ~  0,  cos-1  2u  from  ^  =  0.  However, 
because  of  the  sin#  term  in  Eq.(68),  9  —  0  has  no  contribution  to  the  inte¬ 
gral,  so  we  only  have  to  include  a  stationary  point  u  =  cos  9  —  2v.  Taylor 
expansion  of  around  9  —  cos-1  2 u  with  u  =  cos  9  ~  2u  is  (to  second  order) 
gives 


= 


u  cos  #T 


fu  —  V 


®\u=2v  + 


du 


( u  —  2  v)  + 


i  a2$ 


u— 2v 


2  du 2 


(u  —  2i/)2 


u—2iy 


cos  9X  (69) 


=  COS  9r 


'  r,i  1  ,  9  ^ 

2W+4^75(“-^) 


Then,  Eq.(68)  becomes 


1  —  Au2 

E  —  An2iG(gs^x)Jo  I  A\l - sin  9X 


v 


"°°  -iA  cos  0X  (  ^  3^  , 

e  v  v  4,3/2  (70) 


'2z/ 


And  hnally  we  can  get 


E  =  4t r5/2i1/2G^COB-i2^x)yA'- 


Jo  (jVl  -  4,2  sin6h) 

Vj  COS  ^ 


g  — 2 ?Vy  COS  ^ 


(71) 


-/A  cos  03 


using 


4^/2 


du  — 


2v3/ 


’  2v 


\J iA  cos  0. 


This  equation  shows  that,  if  ~  0,  then  the  wave  propagates  as  E  ~  (£) 
(from  Jq  —  1),  whereas  if  9X  is  sufhciently  large  to  make  the  argument  of  J0 
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large,  then  E  ~  -  ^from  J0  —  y  rsi^  j .  The  transition  occurs  at  ~  - 
from  the  term  inside  of  J0  function  in  Eq.(71).  This  means  that  the  region 
of  weaker  decay,  E  ~  (£)  is  a  cylinder  of  radius  rsinj^  ~  l  about  B 0. 
So,  there  is  a  duct  along  the  B  line  as  in  Fig.  (28).  The  wave  propagates 
intensely  inside  this  duct.  However,  the  word  ”duct”  may  be  misleading, 
because  the  radiating  power  carried  inside  it  is  still  decaying  with  r .  The 
Poynting’s  vector  will  scale  as  E2  ~  £  inside  the  ”duct”,  and  since  its  area 
remains  constant  (7J2),  the  power  it  carries  decays  as  K  Of  course,  outside 

this  duct  the  Poynting’s  vector  decays  as  (^)2,  much  faster.  Now  we  need 

Z 


Figure  28:  A  ”ductw  along  B  line 

to  evaluate  whether  we  can  move  the  J0  term  out  of  the  integral  in  Eq.(67). 
The  J0  term  can  be  moved  out  if  its  argument  y  — 

varies  little  in  an  interval  A u  of  u  about  u  —  2u  such  that  the  exponent  A$f 
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changes  by  about  one  unit.  From  Eq.(69), 


Then, 


Now, 


AA<F  =  A  cos  9X—  — -j-Au2  —  1 

4  u6!1 


A  u  — 


dy_ 

du 


a/AcOsW 
u2  —  2 uu  +  1  A  . 


a/1  —  u2  {u  —  i/)3^2  2 


sin  9r 


(72) 


(73) 


(74) 


And,  in  order  to  move  the  Jo  term  out,  the  condition  below  should  be  satis¬ 
fied: 


dy_ 

du 


\[A 


A?/  _  v  -  sin 9x 
u=2u  VCOS  9X 


<  1 


(75) 


Then,  the  condition  for  9X  (assumed  to  be  small)  is 


9X  < 


(76) 


The  two  limitations  on  9X  can  be  visualized  as  in  Fig. (29).  Above  OC 
(red  line),  the  held  decays  as  ~r  (outside  of  the  ”duct”),  whereas  below  OC, 

Above  OABD  (blue  line),  the  T0  function  varies  too  rapidly  and  must  be  kept 
inside  the  integral,  whereas  below  OABD  it  can  be  moved  outside. 

For  the  far  held  (r  >  p),  there  is  therefore  a  signihcang  region  (OABO)  out¬ 
side  the  high-flux  duct  where  the  J0  function  can  be  moved  out  of  the  integral, 
and  where  Eq.(71)  gives  a  good  representation  of  the  angular  dependence  of 
the  held.  Notice  that  this  could  include  conditions  where  the  argument  of  J0 
is  large,  and  therefore  its  oscillatory  behavior  is  visible.  This  then  overlaps 
with  the  region  where  an  asymptotic  expansion  of  J0  is  valid,  which  turns  out 
to  be  equivalent  to  an  expansion  of  the  exponent  to  second  order  in  </>  —  </>5, 
namely,  the  full  stationary  phase  method. 
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9X  [rad] 


Figure  29:  Jo  in-out  and  propagation  ”duct”  limitation  on  0X 
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The  Earth's  Magnetic  Field 


At  a  large  scale,  The  magnetosphere  shields  the 
surface  of  the  Earth  from  the  charged  particles  of 
the  solar  wind.  It  is  compressed  on  the  day  (Sun) 
side  due  to  the  force  of  the  arriving  particles,  and 
extended  on  the  night  side.  (Image  not  to  scale.) 


magnetic  dipole 


rotational  axis 
(geographic) 


On  a  scale  of  a  few  Earth  radii,  the  field 
approximates  that  of  a  Magnetic  Dipole,  slightly 
tilted  with  respect  to  the  Earth's  axis. 

Along  the  magnetic  shell  crossing  the  equator  at  r 


Equator 


B  ~  R Equator  C0S  ^ 


b  =  be( 


Rt 


D 

v Equatorial 

Re  =  6,130km 


-y 


and  then 

a/4-3cos2  X 


cos6  X 


Br  =3.1x10  Tesla ; 


The  L-Shell  parameter  is  R Equatorid  /  RE  for  a  particular  magnetic  shell. 
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The  Van  Allen  Radiation  Belts 


A  fraction  of  the  high  energy 
particles  from  the  solar  wind, 
plus  others  generated  from 
cosmic  rays,  are  trapped  in 
the  magnetic  bottle  formed 
by  the  local  magnetic  tube. 


Outer 

Radiation 

Bell 


Outer 

Radialic 

Belt 


Rotational 


Inner 

Radiation 

Belt 


Inner 

Radiation 

Belt 


Magnetic  J" 

A  uii.  f  1 


The  outer  belt  contains  mainly  electrons  of  0.1-10  MeV  (energetic  ions  have 
too  large  Larmor  radii  at  these  weak  B-fields). 

Te  inner  belt  contains  electrons  of  up  to  1  MeV,  plus  ions  (mainly  protons) 
of  energies  up  to  400  MeV. 

The  so-called  safe  gap  in  between  is  formed  by  natural  precipitation  of  particles 
Through  pitch  angle  diffusion. 

In  addition,  there  is  a  low-energy  plasma  background,  with  much  higher 
density,  but  much  lower  energy.  This  plasma  is  quasi-neutral.  192 


Effects  on  satellites 


The  bombardment  by  energetic  trapped  particles  damages  solar  cells, 
integrated  circuits  and  sensors  on  board  satellites.  Geomagnetic  storms, 
during  which  teh  radiation  intensifies,  often  disrupts  electronics  on  board  spacecraft. 

These  effects  are  intensified  in  modern  miniaturized  devices,  which  may  have 
Charge  levels  comparable  to  that  of  the  bombarding  ion  and  its  trail. 

Missions  that  need  to  reside  for  long  periods  in  the  Belts  need  to 

(a)  Use  special,  defect-free  components  (“Class  S  electronics”)  that  are  very 
expensive  and  hard  to  find,  as  electronic  manufacturers  dislike  their  low 
production  volume. 

(b)  Use  radiation  hardening  measures,  that  increase  weight  and  cost.  A  3mm 
aluminum  shield  still  allows  exposure  to  2500  rem/year  (Imrem  is  the  dose 
due  to  a  medical  X-ray  exam) 

Manned  missions  are  particularly  critical  in  this  respect.  Apollo  crews  passed 
through  teh  Belts  quickly,  but  high  altitude  orbiting  missions  are  essentially 
precluded. 
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High  altitude  nuclear  explosions 


Another  possible  source  of  trapped  radiation  is  that  generated  in  a  nuclear 
blast  above  the  atmosphere.  Neutrons  in  fission  fragments  (or  free  neutrons) 
can  then  beta-decay  to  protons,  ejecting  high-energy  electrons. 

Crazy  as  it  may  seem,  both  the  US  and  the  former  USSR  did  test  high  altitude 
nuclear  devices: 

US  “Starfish”,  July  9,  1962.  A  1 .4  Mton  device  exploded  at  400  km  altitude  (L=1 .2) 
USSR,  October  22,  1962  .  A  sub-Mton  device,  exploded  at  L=2. 

Both  created  intense  new  radiation  belts  near  the  detonation  L-shells,  that  took 
About  one  year  to  decay.  About  173  of  all  orbiting  satellites  were  damaged. 

The  fear  of  another  such  occurrence,  especially  with  many  new  nuclear  nations 
and  many  more  and  more  valuable  satellites,  has  prompted  research  on  possible 
remedies. 
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E*  MeV ) 


STARFISH  decay  times 


L  (earth  radii) 


Figure  2.12:  Decay  lifetime  (days)  for  Starfish  electrons  (E  >  0.2MeV)  as 
determined  by  the  data  from  the  19 63-33 C  satellite  [9]. 


Particle  spectrum  from  Soviet  nuclear  test 


Fi^urt  1.1:  Tripped  aLccron  ipg^jum  at  viri«ii:  L-i-'duns  c  U  •’viriz;  ch  O- 
bober  23.  l'j*0S!  Scvirt  cuHtfar  decocujcioQ  ']■"/ . 
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Dynamics  of  trapped  particles 

Each  trapped  particle  executes  three  types  of  periodic 
motion: 

Fast  Larmor  motion  about  B 

Intermediate  bounce  motion  between  conjugate  points 
Slow  drift  motion  about  the  magnetic  axis. 

Larmor:  r Larmor  =  — ;  coL  =  —  ~  106  rad  /  s  ( electrons )  or  5,000  rad  /  s  (ions) 

coL  m 

Hence  rL  «  lO  5^  (electrons)  or  0.02  5  (ions) 


Bounce: 


V Pamiiei  =  V  cos  cl  «  10  8  m/  s  (relativistic)  ;  L  «  40,000 km 


Hence  tb  = 


L 


V 


0.4s 


par 


Drift: 


^ drift 


V 


8\2 


par 


(10s) 


^  curvature 


106  X  10 


=  103  m!  s  ;  L  «  80,000 


Hence  rdrift  «8xl04s«l  day 
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Sample  background  plasma  parameters 

In  the  middle  of  the  inner  belt,  para  L=1 .5  (some  3,000km  at  the  Equator) 


Characteristics  of  Natural  Plasma  Environment 

Magnetic  Field 

D  =  9.815x10-*  Gauss 

Plasma  temperature  (daytime) 

T  =  0.4  eV 

Plasma  temperature  nighttime 

T  =  0.2  eV 

Plasma  density 

n  =  I04em-3 

Ion  composition 

H+  =  5x10*  cm-3 

He+  =  iO'Vui-3 

Q+  =  l()2ciu-s 

Neutral  density 

n  =  I04em-3 

Neutral  composition 

H  =  104cm-a 

He  =  10“  cm”3 

O  =  103cm_s 

Larmor  Radius 

re  =  15.3cm 
r[  =  6.57xl02cm 

Debye  Length 

Ad  =  4.70eni 

Collisionality  (mean  free  path) 

Ae  =  135km 

A3  =  192km 
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Earth  r*id  i  i  | 


Variations  with  altitude:  density 


Leg  10  Electron  Density  Midday  13  Jur*e  1996  16'>U  West  Magnct-C  Longitude 


GPSnrhUl  p^jrt 
at  tfrlUHkm  altitude- 


tins  +4,"P 


'5 


X3 


1  HlJjflitu  ITqujttfi 


Ge«^  r.iphic  r 
Eqia^Tilr  (  1«-  E  | 
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Ejrth  tadil(L-»helh 


Figure  2.5:  Electron  densities  (m-3)  described  as  a  function  of  distance  from 
the  earth.  Densities  were  calculated  by  the  Global  Plasmasphere  Ionosphere 
Density  (GPID)  model  [4]. 
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Energetic  electron  fluxes 

Electrons  crossing  per  cm2,  per  second.  Divide  by 
velocity  to  get  particle  density.  Ve  a  c  =  3  x  i010  cm Is 


t  E+D9  j 


Figure  2.4]  High  energy  electron  fluxes  and  energies  as  a  function  of  L- value. 
The  data  reflects  Vampola’s  update  of  the  AE-3  model  using  data  from  the 
Combined  Release  and  Radiation  Effects  Satellite  (CRRES)  [3]. 
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ProinM/Cld^-Sec 


Energetic  proton  fluxes 


Again,  divide  by  velocity  to  get  density,  but  ions  not  relativistic. 


Figure  2.3:  High  energy  proton  fluxes  and  energies  as  a  function  of  L- value 
from  the  A  PS  MIN  model  [3]. 


)3 


A  note  on  relativistic  energy  and  velocity 


The  energy  in  the  previous  data  is  only  the  kinetic  energy  of  the  particle,  ie,  it 
excludes  the  rest  energy: 

2  2 

KE  =  E  -  m0c2 ,  with  E  =  J /?2c2  +  m0c4 ,  and  p 2  =  —  C 

1  —  vie 


Solving  for  the  velocity, 


Which  approaches  c  when  /  m^c2  » \ 


in  the  opposite  limit. 


For  electrons,  mQC2  =  0.51  MeV 
For  protons,  m0c 2  =  9  AO  MeV 


Thus,  most  belt  electrons  are  relativistic,  but  most  protons  are  classical. 
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Magnetic  bottling-Schematic 
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Magnetic  bottling-physics 

When  a  charged  particle  rotates  and  slides  along  a  non-uniform  B-line,  two 
quantities  are  conserved: 


a)  Full  kinetic  energy, 

b)  Adiabatic  invariant, 

Solving  for  the  parallel  energy, 


E  =  E±  +  Eu 
Et 


T-I  l  2 

Eu  =  2WVu 


E  -  jL  tB 


Starting  from  a  region  of  low  B,  where  Eu  >  0  the  parallel  velocity  decreases 

as  B  increases,  and  becomes  zero  at  the  point  where  B  =  Bbounce  =  E /  ju 
Define  the  initial  pitch  angle  a0  =  tan_1(-^-) 

B0  Vu° 

Then  B^ounce 

sm  a0 

Particles  with  a  large  initial  pitch  angle  bounce  early  (at  high  B),  while  particles 
with  small  initial  pitch  angle  penetrate  deep  into  the  higher  B  values. 
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Conditions  for  trapping 


Trapped  particles  are  characterized  by  having  a  bounce  point  above  about 
120  km.,  below  which  there  are  enough  collisions  with  the  neutral  atmosphere 
to  neutralize  or  scatter  the  particle. 


From  r  =  REq  cos2  /i  =  lre  cos2  /I  and  RBounce  *  re 


Substituting  into 


B  _  a/4-3cos2  A  _  1 

BEq  cos6  A  sin2  aEq 


we  get  cos  A  Bounce 


we  then  find  the  pitch  angle  ctEq  corresponding  to  a  given  L-shell: 

^ Ea  —  (  -r^n  ,  a  ^  l / 4  ) 


ZA/z(4-3  ny 


This  gives  27.2  deg  for  L=1 .5,  and  8.4  deg  for  L=3.  Particles  with  a  smaller 
equatorial  pitch  angle  are  quickly  lost  to  the  atmosphere  at  the  bounce  point 

(“Loss  Cone”). 


Any  mechanism  that  scatters  a  trapped  particle  into  an  equatorial 
pitch  angle  smaller  than  the  above  limit  will  therefore  remove  that  particle 
from  the  radiation  belts. 
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Mechanisms  for  particle  precipitation 


Several  natural  effects  have  been  identified  that  do  scatter  particles  into  their 
Loss  Cone: 

Coulomb  Scattering:  Although  the  background  density  is  very  low, 

it  is  not  zero,  and  some  electrons  are  scattered  through  Coulomb  collisions 
with  other  bound  and  free  electrons.  This  dominates  at  the  lower  altitudes, 
were  densities  are  higher. 

Scattering  by  Plasmaspheric  Hiss:  This  is  a  low  frequency  (about  1  kHz) 
broadband  incoherent  radiation  apparently  due  to  cyclotron  resonance  of 
newly  injected  electrons. 

Scattering  by  VLF  Whistler  waves  from  lightning:  These  are  somewhat 
higher  frequency  waves  (1  to  20  kHz)  that  enter  the  ionosphere  near  the 
polar  regions  and  then  propagate  along  B  lines  and  resonate  with  the  gyrating 
electrons. 

In  addition,  it  has  been  realized  since  the  1970’s  that  VLF  radiation  from 
high-power  ground  antennas,  at  discrete  frequencies  in  the  10-30  kHz  range, 
can  at  night  excite  magnetospheric  Whistler  waves  as  well.  The  experimental 
evidence  is  still  only  episodic,  though. 
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Some  model  calculations  of  precipitation  lifetime 


From  Abel&Thorne,  1998,  with  data  from  decay  of  STARFISH: 


Figure  9.  Precipitation  lifetime  calculations  for  500 
kcYr  electrons  for  scattering  due  to  Coulomb  collisions 
(C)*  Coulomb  and  plasmaspheric  hiss  (C/II),  Coulomb* 
plasmasphcrie  hiss  and  lightning-generated  whistlers 
(C/H/W),  and  with  all  scattering  mechanisms  included 
(C/'H/W/YLF).  Observed  decay  rates  are  included  for 
comparison. 


Note  that  the  man-made  VLF  radiation  effect  is  already  very  strong 
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Potential  for  human  intervention 


Since  the  number  of  energetic  particles  involved  is  relatively  small,  there  is 

potential  for  artificially  enhancing  the  precipitation  rate. 

The  drift  around  Earth  (in  hours)  should  spead  any  disturbance  in  longitude. 

Also,  radial  particle  diffusion  between  L-shells  is  slow,  so  one  could  in  principle 

“clean  up”  relatively  thin  shells  for  specific  needs. 

Methods  that  have  been  proposed  include: 

Electrostatic  scattering  of  charges  by  orbiting  thin  bodies  (Tethers)  biased  to 
potentials  of  the  same  order  as  the  particles  to  be  scattered.  With  proper 
geometry,  some  of  the  scattered  particles  will  go  into  the  Loss  Cone. 

VLF  ground  antennas  distributed  strategically  so  as  to  pump  Whistler  radiation 
into  the  desired  L-shells.  This  would  systematize  the  already  observed  effects 
of  operational  (mainly  military)  VLF  installations. 

Orbiting  VLF  or  ULF  antennas,  placed  in  the  specific  targeted  orbital  regions 
to  be  cleaned  up.  Since  the  wavelengths  tend  to  be  hundreds  of  meters  to 
many  km,  a  very  attractive  implementation  would  use  tethers  as  antennas. 


Electrostatic  scattering  of  radiation  belt  particles 


We  next  discuss  some  of  the  details  of  the  electrostatic  scattering  method,  using 
as  the  main  source  the  MS  Thesis  of  C.  Zeineh  at  MIT  (2005). 

A  vertical  tether  in  an  equatorial  orbit  is  a  suitable  geometry.  It  is  dynamically 
stable  due  to  a  Gravity  Gradient  restoring  torque,  and  it  cuts  B  lines  almost 
perpendicularly,  so  that  at  least  some  of  the  scattering  events  can  direct  the 
particles  towards  B,  as  desired. 


The  tether  must  be  at  a  potential  similar  to  the  energy  of  the  particles  to  be 
Affected.  For  most  of  the  study  we  assume  a  1MV  bias. 
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Tether  potentials 


In  order  to  have  a  tether  biased  at  some  potential,  some  sort  of  counter-electrode 
is  required.  This  can  be  a  separate  segment  of  the  same  tether,  or  a  separate  tether. 


Thin  thethers  will  collect  current  in  the  Orbital  Motion  Limit  (OML): 


dl  l2qAV 
—  =  qn  - 

dz  V  m 

For  electrons,  m  =  me  =  0.91  xlCT30  kg  while  protons  are  1840  times  heavier 


Thus,  to  have  equal  current  collected  per  unit  length,  AV 


must  be  1840  times 


smaller  for  electrons  (positive  AV  )  than  for  ions  (negative  AF ).  Since  what  must 
be  the  same  is  really  the  total  current,  this  ratio  could  vary  somewhat  by  assigning 
different  lengths  to  the  two  polarities.  The  result  will  be  a  high  negative  potential 
tether,  ballasted  by  a  moderate  positive  potential  tether  (or  segment). 
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Tether  architectures  (a) 


Single  tether,  no  cathode: 


Positive  tether  provides  ballasting, 
automatically  at  low  potential.  Ob 

Problem:  Net  Loretz  force 
will  drift. 


Positive  Tether  III  Negative  Tether 


Tether  architectures  (b) 


Single  tether,  with  a  cathode: 


No  net  force,  only  a  torque,  can  be 
counterd  by  a  small  angular 
deflection. 

Problem:  No  counter-electrode, 
will  float  at  a  high  positive  potential, 
collect  very  high  current. 


lOkin 


10km 


> 


Tether  architectures  (c) 


Parallel  dual  tethers: 


No  net  force  (some  distortion  only) 
Positive  tether  automatically  at 
low  potential 

Problem:  The  sheath  about  the 
negative  tether  is  of  the  order  of 
200m  radius,  difficult  to  avoid 
interference. 
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Negative  Tether 


Tether  architectures  (d) 


Symmetric  series  tether: 


No  net  force,  only  a  torque  (some  deflection) 


No  interference  due  to  the  sheaths  on  the 
negative  (high  potential)  segments 

Ob 

Low  potential  ballasting  by  center  tether 


Longer  length  provides  additional 
mechanical  stability. 


z 


Selected  design  y 


5km 


10km 


5km 
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Tether  deflection 


Due  to  the  distributed  Lorentz  forces,  a  torque  appears  about  the  cm  of  the 
assembly.  This  torque  is  proportional  to  the  current  I  in  the  cable.  Balancing 
against  the  restoring  gg  torque,  the  deflection  angle  is 

3rofi2 

The  mass  m  of  the  tether  is  in  turn  proportional  to  the  power: 

m  =  ocP  =  al<& 

Substituting,  the  current  cancels  out,  so  the  deflection  ends  up  being 
independent  of  tether  thickness.  In  addition,  B  and  Q2  are  both 

inversely  proportional  to  R 3  so  the  deflection  is  also  independent  of  altitude. 

Using  a  1  MV  bias  and  a  =  20  kg/ kW  we  calculate  a  very  small  deflection: 

6  =  0.02  deg 
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Tether  diameter 

The  current  collected  is  a  necessary  evil,  but  it  can  be  minimized  by  reducing 
the  tether  diameter.  Assume  we  place  an  upper  limit  PMAX  on  the  power 

(product  of  tether  current  and  voltage).  Using  the  OML  formula  for  current,  then 


ma\ 


m 

f 


3  I 

max  | 

-  d 


P  I 

max  j 


\  2e<PT+  eniXL<Z>r  \-2eOr~  rt^L  ^-2[eOv~  f 


rT  <  90.3-  x 


W 


and  for  a  plasma  density  ne  =10* cm  3  =  1 01U m  \  PMAX  =  1  OOkW  =  \0:>W, 


10  „  -3 


MAX 


we  find  a  tether  radius  of  0.9  mm,  or  a  diameter  of  1.8  mm.  Even  accounting 
for  some  extra  current  over  and  above  OML,  this  appears  manageable. 
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Ion  sputtering 


There  is  some  concern  about  the  very  energetic  ions  that  would  bombard  the 
negative  tether.  At  energies  of  a  few  hundred  to  a  few  thousand  eV,  the 
sputtering  yield  (atoms  removed  from  the  metal  per  incident  ion)  can  be 
of  order  unity. 

However,  and  somewhat  paradogically,  the  yield  becomes 

very  small  at  energies  above  lOOkeV.  The  lifetime  of  a  1mm  tether  due  to 

sputtering  turns  out  to  be  many  centuries. 
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Size  of  the  sheath 


The  very  high  negative  potential  about  the  end  tethers  will  produce  a  region 
devoid  of  electrons,  and  where  the  background  plasma  is  therefore  non-neutral 
(a  sheath).  A  smaller  region  devoid  of  ions  will  also  appear  about  the  posive 
middle  tether.The  radius  y  of  the  sheath  can  be  calculated  from  the  implicit 
formula  5 


eO x 

Tr" 


(  Y  V 325 

: 


that  comes  from  numerical  simulations  of  Choiniere,  verified  analytically  by 


Sanmartin  (2007).  Here,  r 


shielding  length, 


is  the  tether  radius,  and  XD  is  de  Debye 


Results  are  shown  below: 


Altitude  (km) 

Plasma 

Debye  Length 

(m) 

Sheath  Radius 
Negative  Tether 

(m) 

Sheath  Radius 
Positive  Tether 

(m) 

2000 

0.040275 

242.68 

1.2346 

4000 

3900 

0.043134 

234.80 

1.1959 

6000 

5100 

0.049326 

220.11 

1.1237 

sooo 

6200 

0.054386 

210.00 

1.0739 

Scattering  calculations 


A  given  high  negative  potential  tether  can  scatter  both  electrons  (repelled  species) 
and  ions  (attracted  species).  An  example  for  repulsion  is  shown  below: 


For  attraction,  the  trajectory  would  instead  “bend  around”  the  tether,  but  the 
deflection  would  be  analogous. 
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Scattering  (continued) 


Some  trajectories  that  were  originally  outside  the  loss  cone,  and  hence  “trapped”, 
will  after  deflection  enter  the  loss  cone, and  the  particle  will  be  counted  as 
instantaneously  lost.  The  loss  cone  in  velocity  space  is  shown  below: 
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Classical  scattering  theory  (in  2D) 


By  using  conservation  of  energy  and  angular  momentum,  classical  theory 
predicts  a  deflection  (given  the  incoming  miss  distance  b  and  velocity  g)  of 


X  -  71  <Pl  Ap=  =  j  — 


dr 


\ 


1  b‘‘ 
l--. 


2q0(r) 


r  mg 

where  <fr(r)  is  the  potential  of  the  scattering  center  (the  tether  in  cross-section) 


and  rm  is  the  radius  of  closest  approach,  which  is  the  solution  of  the  eguation 

r  m?~ 

In  our  case,  we  neglect  the  potential  outside  the  sheath  ( r  >  ),  anc|  the  effect 

of  those  ions  that  are  in  transit  inside  the  sheath.  The  potential  is  then  given  by 
the  classical  electrostatic  solution  between  two  concentric  cylinders: 

=  ,  with  0  =  0  for  r>rm 

In (rx  /  rT ) 


223 


Scattering  (continued) 


The  deflection  integral  must  then  be  split  in  two: 


which  can  be  put  into  the  more  compact  form 


with 
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Some  simple  limits 


Hard  body:  An  electron  with  an  energy  much  less  than  the  tether  potential 
(specifically  with  X»\  almost  “bounces  off’  the  edge  of  the  sheath.  This 
gives  the  simple  limiting  deflection 

-i,  b  \ 

X  —  cos  (— ) 


Soft  body:  In  the  opposite  limit,  a  particle  of  either  polarity  with  very  high  energy 
can  penetrate  deep  into  the  sheath  and  suffer  a  small  deflection  only.  Using  the 
condition  X  «  1  we  now  obtain 


x  =  A  cos-1  >1,_ 


2  q 

mg 2  ln(r,  />>) 


) 


For  all  other  cases,  numerical  integration  is  necessary. 
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Selection  of  the  impact  parameter  b 


First  of  all,  no  incoming  trajectories  can  lie  within  the  loss  cone: 


0  <  6  <  sin 


-i 


v  cos  a. 


Also,  not  every  deflection  places  the  outgoing  trajectory  inside  one  of  the  two 
loss  cones  (along  +Y  or  -Y).  The  condition  on  %  is  y, 


-i 


v 

—  CQStf 


\ 


sin 

u 

_  j-  for  b  >  0 
\+  for  b  <  0 


—  \%  +  -  'T  - 


-] 


v 

—  COStf, 

I . 

Ks 


\ 


/ 


V  =  V  g~  +  V: 


And  the  values  of  b  that  satisfy  it  must  be  found 
from  the  scattering  calculation.  They  are  found  to  lie 
in  the  ranges  (b1,b2)  and  (b3,b4),  and  we  define 
the  total  range 


=  -&2|  +  |i3  -&4 
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The  total  removal  rate 


We  now  have  all  the  elements  to  construct  the  desired  removal  rate.  For  a  tether 
length  lt  ,  and  using  a  “Hollow  Cone”  distribution  (a  Maxwellian  with  the 
particles  within  the  loss  cone  removed),  we  have  for  the  number  of  particles 
removed  per  second: 


.  TCD SON 

™  =  s  ■  J7  J "  J“  '  f  V  ■  fLC  ( v, ,  §,  9)  •  w,  (vr,g,8)  LT  d 6dgdv2 


The  calculations  are  done  numerically  by  discretizing  the  ranges  of  integration. 
C.  Zeineh's  thesis*  explais  the  details  and  gives  interesting  intermediate  results. 

For  reporting,  these  rates  are  divided  by  the  number  density  of  the  energetic 
particles,  yielding  a  quantity  Y  with  units  of  m3 sY 


Applications  of  an  Electrostatic  High-Voltage  Tether  to  Radiation  Belt  Remediation 
M.S.  Thesis,  Depatment  of  Aeronautics/Astronautics,  MIT.  September  2005  227 


Results  for  a  1MV  tether  (a) 


One  first  quantity  of  interest  is  the  the  influx  of  particles  (of  either  polarity)  into  the 
sheaths  of  the  positive  or  the  negative  tether. 


Altitude  (Ian) 

Volume  Influx  (mV1) 
Positive  Tether 

2000 

i— 1 

© 

T*-l 

* 

FO 

i 

un 

¥=l 

2.9786*  i0u 

4000 

1.4678*  i0u 

6000 

1.3792*  i0u 

2.7016*  i0u 

8000 

1.3181*  iO11 

2.5775*  i0lj 

Of  course,  the  main  interest  is  in  the  negative  tether,  which  has  the  high  potential 
and  therefore  the  higher  flux. 
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Results  for  a  1MV  tether  (b) 


Next  we  report  the  volumetric  scattered  fluxes,  and  the  ratio  of  scattered  to 
incoming  flux,  or  “scattering  efficiency”: 


Altitude  (km) 

Volume  Scattering  Flux 
Positive  Tether  (mV1) 

Volume  Scattering  Flux 

■J  4 

Negative  Tether  (m  s  ) 

2000 

8,6978*  lO9 

1.8130*  iO12 

4000 

3.9829*  *0y 

8.6175*  iO1' 

6000 

1.3792*  iO11 

4.4095*  1011 

8000 

1.3181*  40u 

1.5059*  i0w 

Altitude  (km) 

Efficiency  (%) 

Positive  Tether 

Efficiency  (%) 

Negative  Tether 
— 

2000 

5.74 

6.09 

4000 

2.71 

2.99 

6000 

1.45 

1.63 

8000 

0.02 

0.06 

The  higher  efficiencies  at  the  lower  altitudes  reflect  the  fact  that  the  loss  cone 
is  wider  at  these  altitudes.  22 


“Remediation”  time 


The  rate  of  decrease  of  particle  density  being  proportional  to  the  density  itself, 
the  evolution  of  this  density  (or  of  the  associated  flux)  turns  out  to  be  exponential 

_L.t 

n  =  n0e  v 


where  V  is  the  volume  of  the  region  being  “remediated”.  The  time  to  achieve 
a  particular  decrease  ratio  follows  as 


{ 


\ 


y 

i 


ffwal  J 


The  volume  V  is  very  roughly  estimated  as  that  comprised  between  two  toroids 
whose  major  and  minor  radii  are  both  half  the  equatorial  radius  of  the  particular 


L-shells  that  bound  the  “remediation  region”: 


=  \x2  ~  R 

4 


3 


eq. inner 


) 


And  further,  we  assume  the  “remediation  region”  is  bounded  by  the  orbital  radius 

of  the  tether  plus/minus  a  distance  D,  that  we  equate  with  the  10km  teher  length. 
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Remediation  time  results 


Setting  as  the  target  a  reduction  to  1/10  of  the  original  energetic  particle  flux, 
we  obtain  the  following  results  : 


Altitude 

T  (sec) 

T  (years) 

2000 

1.3136*  i0,; 

0.4166 

4000 

6.0330*  i07 

1.9130 

6000 

2.0643*  10® 

6.6458 

8000 

1.5722*  1010 

498.54 

These  results  assume  no  diffusion  from  adjacent  L  shells.  To  see  the  possible  effect 
of  this  diffusion,  suppose  we  increase  the  distance  D  to  100  km  above  and  below 
the  ends  of  the  tether.  We  no  calculate  the  longer  times 


Altitude 

T  (sec) 

T  (years) 

2000 

1.3135*  30s 

4.1652 

4000 

6.0331*  108 

19.131 

6000 

2.0643*  10y 

65.459 

8000 

1.5722*  i0u 

4985.5 

The  mission  appears  feasible  only  for  the  lower  L-shells,  to  about  L=1.^.5 


Conclusions  for  the  electrostatic  scattering  method 


-  No  apparent  technology  barrier 

-  Can  work  for  either  electrons  or  ions 

-  Restricted  to  the  inner  belts,  probably  to  its  lower  edges 

-Fairly  high  power  required  to  maintain  the  tethers  at  high  potential 
(of  the  order  ot  100  kW  for  1  mm  diameter) 

-Parametric  studies  needed  for  other  geometries,  multiple  tethers,  etc. 


232 


Methods  based  on  induced  pitch  diffusion 


The  eletrostatic  scattering  method  relies  on  deflecting  a  fraction  of  the  affected 
particles  into  the  loss  cone  all  in  one  very  energetic  encounter  with  the  charged 
tether. 

The  difusion  methods,  to  be  considered  next,  rely  on  a  multiplicity  of  “gentle” 
deflections,  that  accumulate  as  the  particle  travels  repatedly  back  and  forth 
along  its  guiding  magnetic  line  (actually  surface,  considering  the  drift). 

These  methods  can  work  because  the  scattering  agent  is  a  wave  field  that  can 
be  made  to  essentially  “fill”  the  same  magnetic  shell  as  the  particles.  However, 
the  resonance  conditions  needed  for  efective  interaction  occur  only  over  a  few 
percent  of  the  particle's  path. 

The  physics  of  this  interaction  is  fairly  complex  in  detail,  and  has  been  studied 
by  space  physicists  over  the  past  40  or  so  years  in  order  to  explain  a  multitude 
of  magnetospheric  and  ionospheric  effects.  Only  recently  however  has  the 
artificial  excitation  of  these  effects  been  seriously  considered. 
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Wave  fields 


A  wave  propagating  through  the  background  plasma  will  carry  oscillating  electic 
and  magnetic  fields,  Ew  and  Bw.  Normally  Bw  «  B0  the  background  B  field. 

These  fields  are  related  to  each  other  through  one  of  Maxwell's  equations: 

as 


dt 


=  -V  x  E, 


For  a  planar  wave,  every  perturbation  quantity  a  will  vary  in  the  Fourier  form 

a(r,t )  =  ${{a(k,co)exp(i(k.r  -  cot)]} 

where  the  frequency  CO  is  a  certain  function  of  the  components  of  the  wave 

vector  k  as  determined  by  the  properties  of  the  plasma.  The  phase  (the  exponent) 
will  be  constant  for  an  observer  moving  at  the  phase  velocity, 

£  V 

In  Fourier  components,  we  have  a)B  =  kxE  which  shows 


(a)  E  and  B  are  perpendicular  to  each  other,  and 

E  _  co  _ 

B~  k~V 


(b)  In  magnitude, 


* 
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Wave-particle  interaction 

The  wave  fields  exert  a  force  on  a  particle  of  velocity  v  and  charge  q 

F  =  <l(Ew  +  vxBw) 

-  Since  the  magnetic  force  is  perpendicular  to  the  velocity  it  cannot  add  energy 

to  it,  only  deflect  its  trajectory. 

-  The  ratio  of  the  magnetic  and  electric  forces  is  Ll  =  sm g  =  _Lsin(~ 

FE  Ew  ^ 

where  <P  js  the  angle  of  the  velocity  and  the  wave  B  field. 

-  The  velocity  ratio  —  =  n  is  called  the  index  of  refraction  of  the  plasma. 

For  the  whistler  waves  of  interest,  n  is  of  the  order  of  20  or  more,  meaning 
these  waves  are  much  slower  than  the  speed  of  light  c  in  vacuum.  On  the  othei 
hand,  the  particle  velocity  v  is,  at  least  for  electrons,  close  to  c,  as  we  saw. 
Hence 

(a)  The  main  wave-particle  interaction  is  magnetic  deflection 

(b)  The  energy  of  the  particles  is  affected  only  weakly  by  tese  waves 


235 


The  need  for  a  resonance 


-  The  basic  magnetic  field  B0  changes  the  velocity  direction  appreciably  in  a  time  of 
the  order  of  1/Qe  ;  Cle  =  eB0  /  m  The  wave  field,  being  weaker  by  the  ratio 

Bw/  B0&  10“  will  require  a  time  proportionally  longer  for  a  similar  deflection. 

-  When  the  frequencies  of  the  wave  and  of  the  gyrating  particle  are  unrelated,  the 

magnetic  forces  will  change  direction  rapidly  and  more  or  less  at  random,  with 
not  much  net  effect.  But  when  they  coincide,  or  the  wave  frequency  is  a  multiple 
of  the  gyration  frequency,  (a  resonant  condition)  the  effect  can  accumulate  for 
a  reasonably  long  time,  and  lead  to  measurable  deflections.  In  addition, 

(a)  The  wave  frequency  must  be  Doppler-shifted,  so  that  it  is  the  frequency 
seen  by  the  advancing  particle. 

(b)  The  gyration  time  must  be  modified  by  the  relativistic  dilation. 

(c)  The  wave  fields  must  rotate  in  the  same  sense  as  the  particle 

-  The  resonance  condition  is  then:  co  -  Lv»  = - - 

I -  7 

with  r=0,  -1 ,  -2....,  and  l//  =  ^l-v2/c2  As  noted  the  wave  frequency 
must  be  as  corresponds  to  the  parallel  wavevector  ku 
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The  Pitch  Diffusion  approximation 

Even  with  the  help  of  the  resonance,  the  deflections  will  happen  at  random, 
depending  on  the  phase  difference  at  the  start  of  the  resonant  interaction.  As 
the  particle  bounces  back  and  forth  between  conjugate  points,  these  random 
deflections  will  produce  a  diffusion  of  the  particle  population  in  the  pitch  variable. 


Following  these  ideas,  theorists  have  constructed  the  formulas  that  determine  the 
diffusivity  in  pitch.  These  are  different  for  each  of  the  resonances  possible 
(r=0,-1,-2,  ...).  They  involve  averaging  along  the  guiding  magnetic  lines,  as  well 
as  over  each  Larmor  cycle,  and  they  depend  on  the  specific  equatorial  pitch,  a0 


As  an  example,  for  the  r<0  resonances,  the  pitch  diffusivity  (in  rad 2  Is  )  is 


fnr  ' 

(- 

A  f 

\  TC 

Bz 

w 

P\\,m 

lao  J 

r(a0 

)cos2  a0  ^ 

1  2 

B2 

Y 

p\\ 

W„  cos  a  cos7  X  dX 


The  details  are  explained  for  instance  by  J.  M.  Zorrilla  *  in  his  report. 


(*)  J.  Manuel  Zorrilla  Matilla,  “Radiation  Belt  Remediation  Methods”.  Graduating 
Report  for  Supaero,  (done  at  MIT,  Spring  2006).  ^37 


The  pitch  diffusion  equation 

For  a  given  particle  energy  E  and  magnetic  shell,  we  can  define  the  distribution 
function  in  pitch,  f0(E,L,t, oc0)  as  the  number  of  particles  per  unit  pitch  width 

and  unit  volume,  at  time  t  for  which  the  equatorial  pitch  is  a0 


The  distribution  function  is  governed  by  the  time-dependent  pitch  diffusion  equation 


dt 


sin2aor(a0 


j 


V  V 


The  boundary  conditions  are  (a)  Symmetry  about  zero  pitch,  and  (b)  Zero 
population  at  the  edge  of  the  hollow  loss  cone.  Since  these  are  homogeneous 
conditions  on  a  linear  PDE,  any  non-trivial  solution  must  be  an  eigenfunction 
of  the  problem. 

For  solution,  separation  of  variables  is  used:  /o(^ao)  =  ^k(ao) 

and  the  separation  constant  \/  ?  plays  the  role  of  an  eigenvalue: 


1  dNjt) 
N(t)  dt 


1  d 

sin2aor(a0  )g(a0 )  daQ 


f 

sin 

V 


2a0r(a0)£>ao(a0) 


dgio-S 

<ia0  j 


T 
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Numerical  solution 


Several  authors  have  solved  this  separated  equation  by  iteration.  Zorrilla  chose 
to  discretize  the  variables  and  convert  the  formulation  to  a  matrix  eigenvalue 
problem.  Some  care  had  to  be  exercised  to  select  the  lowest  eigenvalue, 
corresponding  to  the  fastest  depletion  time. 


239 


